Frontiers  in  Anisotropic  Shock- Wave  Modeling 

by  Alexander  A.  Lukyanov  and  Steven  B.  Segletes 


ARL-TR-5878 


February  2012 


Approved  for  public  release;  distribution  is  unlimited. 


NOTICES 


Disclaimers 


The  findings  in  this  report  are  not  to  be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated 
by  other  authorized  documents. 

Citation  of  manufacturer’s  or  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  thereof. 
Destroy  this  report  when  it  is  no  longer  needed.  Do  not  return  it  to  the  originator. 


Army  Research  Laboratory 

Aberdeen  Proving  Ground,  MD  21005-5066 

ARL-TR-5878  February  2012 


Frontiers  in  Anisotropic  Shock- Wave  Modeling 

Alexander  A.  Lukyanov 

Abingdon  Technology  Centre,  Schlumberger,  Abingdon,  0X14  1UJ,  UK 

Steven  B.  Segletes 

Weapons  and  Materials  Research  Directorate,  ARL 


Approved  for  public  release;  distribution  is  unlimited. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering 
and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  the  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to 
comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

13.  SUPPLEMENTARY  NOTES 

authors’  email:  aaluk@mail.ru,  steven.b.segletes.civ@mail.mil 

A.  Lukyanov  contact  info:  Tel.:  +44  07840355383,  Fax:  +44(0)  1234  758217. 

14.  ABSTRACT 

Studies  of  anisotropic  materials  and  the  discovery  of  various  novel  and  unexpected  phenomena  under  shock  loading  has 
contributed  significantly  to  our  understanding  of  the  behavior  of  condensed  matter.  The  variety  of  experimental  studies  for 
isotropic  materials  displays  systematic  patterns,  giving  basic  insights  into  the  underlying  physics  of  anisotropic  shock-wave 
modeling.  There  are  many  similarities  and  significant  differences  in  the  phenomena  observed  for  isotropic  and  anisotropic 
materials  under  shock-wave  loading.  Despite  this,  the  anisotropic  constitutive  equations  must  represent,  mathematically  and 
physically,  the  generalization  of  the  conventional  constitutive  equations  for  isotropic  material  and  reduce  to  the  conventional 
constitutive  equations  in  the  limit  of  isotropy.  This  report  presents  the  current  state  of  the  art  in  the  experimental  and  theoretical 
developments  of  this  fascinating  field. 


15.  SUBJECT  TERMS 

anisotropic  material,  anisotropic  plasticity,  shock  waves,  equation  of  state,  stress  decomposition 

17.  LIMITATION  18.  NUMBER 
OF  ABSTRACT  OF  PAGES 

UU  72 


19a.  NAME  OF  RESPONSIBLE  PERSON 

Steven  B.  Segletes 

19b.  TELEPHONE  NUMBER  ( Include  area  code) 

410-278-6010 


16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

Unclassified 

Unclassified 

Unclassified 

1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

February  2012  Final 

4.  TITLE  AND  SUBTITLE 

Frontiers  in  Anisotropic  Shock- Wave  Modeling 


6.  AUTHOR(S) 

Alexander  A.  Lukyanov 
Steven  B.  Segletes 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Abingdon  Technology  Center,  Schlumberger,  Abingdon,  OK14  1UJ,  UK  ; 
U.S.  Army  Research  Laboratory,  ATTN:  RDRL-WMP-C,  Aberdeen  Proving 
Ground,  MD  21005-5066 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


3.  DATES  COVERED  (From  -  To) 
January  201 1-October  201 1 
5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


5d.  PROJECT  NUMBER 

AH80 

5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ARL-TR-5878 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 


11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


11 


Standard  Form  298  (Rev.  8/98) 
Prescribed  by  ANSI  Std.  Z39.18 


Contents 


List  of  Figures  v 

List  of  Tables  vii 

Acknowledgment  viii 

1.  Introduction  1 

1.1  Experimental  Framework .  2 

1.2  Linear  Anisotropic  Constitutive  Relationship .  5 

2.  Isotropic  Shock- Wave  Modeling  7 

3.  Anisotropic  Shock- Wave  Modeling  12 

3.1  Characteristics-based  Solver .  12 

3.2  Nonlinear  Anisotropic  Elastic  Incremental  Formalism .  16 

3.3  Conventional  Decomposition  of  Anisotropic  Equations .  18 

3.4  Deviatoric  Constitutive  Relationship  for  Anisotropic  Materials . 20 

3.5  Generalized  Decomposition  of  Anisotropic  Constitutive  Relationship  Suitable  for 

Shock- Wave  Modeling . 27 

4.  Anisotropic  Plasticity  Flow  33 

4.1  Associated  and  Non-associated  Anisotropic  Plasticity . 33 


iii 


4.2 


Dislocation-based  Plasticity 


37 


5.  Early  Anisotropic  Implementations  in  Hydrocodes 


38 


6.  Conclusion 

40 

7.  References 

42 

Distribution  List 

57 

IV 


List  of  Figures 


Figure  1.  Schematic  diagram  of  the  experimental  target  assembly .  6 

Figure  2.  Cold-compression  and  shock-Hugoniot  curves  for  aluminum  (a)  to  2.4  and  (b) 

11  megabars.  Note  that  cold-compression  data  are  filled  symbols  and  Hugoniot  data 
are  open  symbols .  11 

Figure  3.  Measured  (solid  lines)  and  simulated  (dashed  lines)  longitudinal  stress  histo¬ 
ries  for  LiF  single  crystals  shocked  along  the  (100)  orientation.  The  experimental  data 
(Asay  et  al.)(73)  and  simulations  (Winey  and  Gupta)  (29)  are  for  a  6061-T6  aluminum 
plate  impacting  the  LiF  crystals  at  a  velocity  of  ca.  0.34  km/s,  where  the  target  crystals 
are  backed  by  a  quartz  stress  gauge.  The  stress  is  measured  at  the  interface  between  the 
LiF  crystal  and  the  quartz  gauge.  The  numbers  above  the  curves  indicate  the  sample 
thickness  in  mm.  Time  is  relative  to  the  moment  of  impact.  Reprinted  with  permission 
from  Winey,  J.M.,  Gupta,  Y.  M.,  J.  Appl.  Phys.  99,  023510,  (2006).  Copyright  2006, 
American  Institute  of  Physics .  19 

Figure  4.  The  Kevlar/Epoxy  IFPT  simulated  and  experimental  back  surface  velocities  for 
572,  788,  and  1015  m/s.  The  experimental  data  Kevlar/Epoxy  materials  recovered  after 
flyer  plate  testing  were  taken  from  C.  J.  Hayhurst,  Multi-Physics  Analysis  of  Hyperve¬ 
locity  Impact:  Successes  and  Challenges ,  FENET  Presentation  (Noordwjik,  2003)  (133).  26 

Figure  5.  Representative  experimental  gauge  traces  from  the  through-thickness  orienta¬ 
tion  at  the  0  mm  position  and  at  the  back  surface,  respectively  (see  Millett  et  al.)(15). 

The  specimen  was  3.8  mm  thick.  The  impact  conditions  were  a  5  mm  dural  flyer  at 
V  =  504  m/s.  The  dotted  curve  is  the  numerical  data  obtained  using  proposed  damage 
model;  the  solid  curve  is  the  experimental  data . 30 

Figure  6.  Experimental  data  Ug-up  for  the  carbon-fiber-composite  material,  showing  the 
variation  with  specimen  thickness  (experimental  data  obtained  by  Millett  et  al.)(15). 

The  dotted  curve  is  calculated  using  experimental  data  for  U[)  -up;  the  solid  curve  is 
calculated  using  numerical  simulation  based  on  the  material  model  equations  66-73 .  31 


v 


Figure  7.  Representative  experimental  gauge  traces  from  the  fiber  0°  orientation  (see  Mil- 
lett  et  al.)(15).  The  specimen  was  10  mm  thick.  The  impact  conditions  were  a  5  mm 
copper  flyer  at  V  —  936  m/s.  According  to  Hereil  et  al.  (11),  the  precursor  was  due  to 
a  high  velocity  wave  transmitted  along  the  fibers  orientated  in  the  shock  axis,  while  the 
main  shock  was  transmitted  through  the  matrix . 32 

Figure  8.  Back-surface  gauge  stress  traces  (longitudinal  direction)  from  plate-impact  ex¬ 
periments  vs.  numerical  simulation  of  stress  (PMMA)  waves  for  plate  impact  tests  (im¬ 
pact  velocities  450  and  895  m/s)  -  target  AA7010  T6 . 36 

Figure  9.  Back-surface  gauge  stress  traces  (transverse  direction)  from  plate-impact  experi¬ 
ments  vs.  numerical  simulation  of  stress  (PMMA)  waves  for  plate  impact  tests  (impact 
velocities  450  and  895  m/s)  -  target  AA7010  T6 . 36 


vi 


List  of  Tables 


Table  1.  EOS  parameters  for  select  materials. 


Acknowledgment 


The  author  (Lukyanov)  thanks  Prof.  E.  Romenski  for  many  useful  suggestions  regarding 
shock-wave  modeling.  The  discussions  regarding  the  shock-wave  experiments  on  a  composite 
with  Dr.  J.C.F.  Millett  during  the  meetings  at  Cranfield  University  are  also  greatly  appreciated. 
Both  authors  would  like  to  thank  Dr.  Bryan  Love  for  his  thorough  and  thoughtful  review  of  this 
report. 


1.  Introduction 


Investigation  of  anisotropic  material  behavior  (such  as  aluminum  alloys,  single  crystals, 
composite  materials,  and  laminated  plates)  has  found  significant  interest  in  the  research 
community  due  to  the  widespread  application  of  anisotropic  materials  in  aerospace  and  civil 
engineering  problems.  For  example,  aluminum  alloys  are  one  of  the  main  materials  in  the 
construction  of  modem  aircraft  and  rockets.  The  strain-rate-dependent  mechanical  behavior  of 
anisotropic  material  (e.g.,  aluminum  alloys)  in  air  and  space  vehicles  is  important  for  applications 
involving  impact.  These  applications  cover  a  wide  range  of  situations  such  as  crashworthiness 
and  protective  armours  in  air  and  space  vehicles  and  other  applications.  Since  shock-wave  and 
high-strain-rate  phenomena  are  involved  in  many  physical  phenomena,  we  are  interested  in 
understanding  the  material  mechanical  properties  under  these  non-trivial  conditions.  This  report 
presents  the  current  state  of  the  art  in  the  experimental  and  theoretical  developments  of  the 
shock-wave  propagation  in  anisotropic  solids,  and  more  specifically,  for  anisotropic  elastic-plastic 
solids. 

This  is  a  subject  that  has  received  considerable  attention  in  the  isotropic  solid-state  physics  and 
mechanics  literature  in  recent  decades  (e.g.,  Wackerle  [7];  Zel’dovich  and  Raizer  [2];  Davison  and 
Graham  [2] ;  Eliezer  et  al.  [ 4 ] ;  Asay  and  Shahinpoor  [5] ;  Meyers \6\ ;  Drumlieller  [7]).  The  impact 
loading  of  anisotropic  materials  is  a  subject  that  has  received  considerable  attention  for  both  low 
and  high  velocity  in  the  recent  monographs  on  composite  materials  (e.g.,  Abrate  [8],  Reid  and 
Zhou  [9],  Ryan  et  al.  [70])  and  in  recent  research  papers  on  composite  materials  (e.g.,  Hereil  et 
al.  [11],  Bordzilovsky  et  al.  [12],  Dandekar  et  al.  [13],  Espinosa  et  al.  [14],  Millett  et  al.  [75], 

Akta§  et  al.  [16],  Hazell  and  Appleby-Thomas  [77],  Garcfa-Castillo  et  al.  [78],  Tekalur  et  al.  [79]), 
research  papers  on  metals  (e.g..  Butcher  [20],  Johnson  and  Barker  [27],  Stevens  and  Tuler  [22], 
Rosenberg  et  al.  [23,  24],  Gray  III  et  al.  [25,  26],  Rubin  [27]),  and  research  papers  on  single 
crystals  (e.g.,  Winey  and  Gupta  [28-30]). 

To  describe  the  anisotropic  material  response  under  shock  loading,  the  following  general  aspects 
need  to  be  investigated:  appropriate  constitutive  equations  to  describe  the  strength  effect  and 
equations  of  state  (EOSs)  to  describe  the  hydrodynamic  behavior.  Mechanical  yielding  and 
strength  behavior  in  shock  waves  show  complexities  that  are  not  understood  yet,  especially  in 
anisotropic  materials.  Some  approaches  have  been  made  by  Johnson  et  al.  (31-33),  Segletes  (34), 
O’Donoghue  et  al.  (25),  Anderson  et  al.  (36),  Winey  and  Gupta  (28-30),  and  Lukyanov  (2 7-41)  to 
describe  the  constitutive  relationships  for  an  anisotropic  material.  The  determination  of  the  EOS 
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for  anisotropic  materials  is  an  important  problem  in  metallurgy,  geophysics,  aerospace,  and  also 
in  other  areas  where  the  behavior  of  anisotropic  materials  at  high  pressures  is  of  interest. 

1.1  Experimental  Framework 

Modem,  high-resolution  methods  for  monitoring  the  stress  and  particle  velocity  histories  in  shock 
waves  and  equipment  have  been  created  (e.g.,  Barker  and  Hollenbach [42];  Kanel  [45];  Kanel  et 
al.  [44];  Millett  and  Bourne  [45];  Bourne  and  Stevens  [46];  Bourne  [47];  Gu  and 
Ravichandran  [48]);  numerous  investigations  into  the  mechanical  properties  of  different  classes  of 
materials  have  been  undertaken  (e.g.,  Meyers  [6];  Gu  and  Ravichandran  [48];  Steinberg  [49]; 
Johnson  et  al.  [50];  Kanel  et  al.  [57];  Millett  et  al.  [52];  Lopatnikov  et  al.  [55];  Zaretsky  et  al.  [54]; 
Gebbeken  et  al.  [55];  Bronkhorst  et  al.  [56]),  and  numerous  phenomenological  as  well  as 
microscopic  models  have  been  developed  (e.g.,  Wallace  [57];  Swegle  and  Grady  [58]; 

Steinberg  [49];  Meyers  [6];  Kanel  et  al.  [59];  Nellis  et  al.  [ 60] ;  Bourne  and  Gray  III  [67];  Kruger  et 
al.  [62];  Chijioke  et  al.  [65];  Boidin  et  al.  [64];  Petit  and  Dequiedt  [65]).  However,  in  spite  of  a 
perfectly  adequate  general  understanding,  experimental  methodology,  and  theory,  material 
models  do  not  agree  in  detail,  especially  for  anisotropic  materials. 

In  the  past,  very  few  experimental  results  were  published  in  the  literature,  either  for  anisotropic 
metals  or  composite  materials.  Only  a  few  examples  can  be  presented.  For  example,  Whittier 
and  Peck  (66)  measured  the  response  due  to  low  impact  and  studied  the  dispersion  of  the  wave  in 
unbounded  layered  media,  not  a  laminated  plate.  Tauchert  and  Guzelsu  (67)  and  Asay  et  al.  (68) 
used  ultrasonic  methods  to  study  dispersion  and  measure  the  “effective”  elastic  response  of 
unbounded  composites.  The  transient  strain  histories  measurements  were  performed  by 
Mortimer  et  al.  (69)  for  fabricated  graphite -epoxy  plates  (cross-ply  layup  and  angle  plies). 
Chhabildas  and  Swegle  (70,  71)  performed  pressure-shear  experiments,  where  the  coupled 
longitudinal  and  transverse  motion  generated  by  the  normal  impact  of  Y-cut  quartz  is  transmitted 
into  an  X-cut  quartz,  X-cut  quartz,  and  alumina- filled  epoxy  samples. 

Jones  and  Mote  (72)  presented  experimental  studies  on  shock  propagation  in  a  single  crystal  in 
copper  including  constitutive  equations  describing  elastic  precursor  decay  for  wave  propagation 
in  the  (100),  (110),  and  (111)  directions  in  fee  single  crystals.  Johnson  et  al.  (31)  proposed 
constitutive  equations  describing  elastic  precursor  decay  for  longitudinal  plane- wave  propagation 
in  fee,  bcc,  and  rocksalt  structures  with  the  wave  propagation  in  the  (100),  (110),  and  (111) 
directions.  Calculated  theoretical  results  (Johnson  et  al.  [5i])  are  compared  with  the 
experimental  data  on  precursor  amplitudes  for  single-crystal  copper  (fee),  tungsten  (bcc),  sodium 
chloride  (NaCl)  (rocksalt),  and  lithium  fluoride  (LiF)  (rocksalt).  Asay  et  al.  (73)  performed 
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experimental  studies  on  shock  propagation  along  a  (100)  direction  in  single-crystal  LiF  and 
showed  that  elastic  precursor  decay  is  critically  dependent  on  the  origin  of  the  sample.  Also, 
Gupta  (74)  presented  experimental  data  for  shock  propagation  along  the  (111)  crystallographic 
direction  in  single-crystal  LiF,  with  the  observed  elastic  response  up  to  30  kbar  experimental 
bounds  and  including  shock  propagation  data  along  the  (100),  (110),  and  (111)  directions 
(Gupta  [75]).  The  dynamic  elevated-temperature  elastic  properties  of  single-crystals  and 
poly  crystalline  aluminum  were  studied  by  employing  a  laser  pulse  technique  to  produce 
propagating  stress  pulses  in  slender  rods  (Swearengen  et  al.  [76]).  The  dynamic  elastic-plastic  of 
single  crystal  of  the  secondary  explosive  pentaerythritol  tetranitrate  (PETN)  was  studied  by 
Halleck  and  Wackerle(77)  at  input  shock  strengths  between  0.6  GPa  for  crystals  with  surfaces  cut 
perpendicular  to  the  (110)  and  (001)  directions.  Later,  Dick  and  Ritchie  (78)  presented  the 
measurements  for  the  elastic  precursor  shock  strength  of  pentaerythritol  tetranitrate  explosive 
crystals  with  (100),  (101),  (110),  and  (001)  orientations.  The  measured  precursor  amplitudes 
were  0.38,  0.58,  0.98,  and  1.22  GPa,  respectively,  for  the  four  orientations. 

However,  this  situation  has  changed  recently  for  composite  materials  and  single  crystals,  i.e.,  a 
number  of  monograph  and  research  papers  were  published  (e.g.,  Abrate  [8];  Reid  and  Zhou  [9]; 
Millett  et  al.  [75];  Ryan  et  al.  [70];  Akta§  et  al.  [7(5];  Hazell  and  Appleby-Thomas  [77]; 
Garcfa-Castillo  et  al.  [78];  Tekalur  et  al.  [79];  Nurick  et  al.  [79];  Iqbal  et  al.  [80];  Enfedaque  [87]; 
Panahia  et  al.  [82];  Jackson  and  Shukla  [85]).  The  major  experimental  methods  were  presented 
by  Zhu  and  Lu  (84)  with  a  brief  description  of  the  corresponding  experimental  devices  such  as 
ballistic  pendulums  and  a  wide  range  of  sensors  to  measure  impulse,  pressure,  acceleration,  and 
displacement  of  structures.  The  energy  profile  diagrams  and  associated  load-deflection  curves  for 
a  number  of  composite  materials  (e.g.,  unidirectional  glass/epoxy  laminates,  carbon  fiber 
reinforced  plastic  [CFRP]  laminates,  woven  carbon  fiber/epoxy  laminates)  were  obtained  (Found 
and  Howard  [85];  Delfosse  and  Poursartip  [8(5];  David- Westa  [87];  Akta§  et  al.  [7(5];  Iqbal  et 
al.  [80];  Enfedaque  [87]).  New  experimental  data  have  been  obtained  for  single  crystals 
(e.g.,  Asay  [88]).  A  good  overview  of  new  development  in  the  physical  chemistry  of  shock 
compression  was  presented  by  Dlott  (89). 

For  many  years,  it  has  been  assumed  in  the  shock-wave  community  that  the  response  of  materials 
to  high  intensity  shock  loading  is  isotropic,  and  only  recently  has  anisotropy  in  the  shock 
response  attracted  the  attention  of  researchers  (e.g.,  Gupta  [75],  Dick  and  Ritchie  [78],  Gray  III  et 
al.  [25,  26]).  It  was  shown  in  an  investigation  of  cold  rolled  and  annealed  zirconium  (see  Gray  III 
et  al.  [25])  that  the  value  of  stresses  varies  in  different  directions  in  the  quasi-static  test  and  plate 
impact  test.  Gray  III  et  al.  (26)  showed  that  under  shock  loading  conditions  (one-dimensional 
strain  space),  the  variation  of  the  Hugoniot  elastic  limit  (HEL)  or  the  yield  strength  of  annealed 
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zirconium  was  consistent  with  the  quasi-static  experimental  data.  Previously,  Butcher  (20) 
worked  on  aluminum  alloy  6061-T6  and  predicted  that  spall  strength  should  vary  in  accordance 
with  the  one-dimensional  yield  strength  depending  on  material  orientation.  Chhabildas  and 
Swcglc  (70)  presented  a  technique  to  detect  a  0.2  GPa  shear  wave  in  6061-T6  aluminum  alloy  at 
0.7  GPa  longitudinal  stress.  Gupta  (75)  presented  shock  propagation  data  along  the  (100),  (110), 
and  (111)  directions  in  LiF  crystals  where  a  marked  anisotropy  in  wave  profiles  and  in  dynamic 
compressive  strengths  is  observed.  Here  it  is  important  to  mention  the  work  done  by  Johnson  and 
Barker (20),  Stevens  and  Tuler(22),  and  Rubin (27)  on  aluminum  alloy  6061-T6  and  Rosenberg  et 
al.  (23,  24)  on  the  shock  response  of  the  alloy  2024-T86. 

By  their  very  nature,  two-dimensional  fiber  composites  are  highly  anisotropic.  Eden  et  al.  (90) 
used  high-speed  photography  to  investigate  a  quartz-phenolic  resin  composite.  Bordzilovsky  et 
al.  (12)  examined  the  effects  of  orientation  to  the  shock  axis  of  a  unidirectional  aramid 
fiber-epoxy  composite,  with  the  orientation  of  the  fibers  ranging  from  5°  to  90°  to  the  shock  axis. 
Work  by  Hereil  et  al.  (11)  on  a  three-dimensional  carbon-carbon  composite  (in  this  case,  the  fibers 
and  binder  are  orientated  orthogonally)  were  in  agreement  with  previous  studies  of  Eden  et 
al.  (90).  Millett  et  al.  (15)  studied  the  effect  of  fiber  orientation  on  the  shock  response  of  a 
two-dimensional  carbon  fiber-epoxy  composite  (through  thickness  [fibers  normal  to  the  impact 
axis]  and  fiber  at  0°  [fibers  parallel  to  the  impact  axis])  using  the  technique  of  plate  impact. 

A  technique  for  the  study  of  the  behavior  of  metals  under  extreme  conditions  is  the  planar  plate 
impact  test  (one-dimensional  plane-strain  shock-wave  propagation).  To  describe  the  dynamic 
response  of  different  materials  under  shock  loading,  the  methodology  based  on  the  plate  impact 
test  is  fundamental  in  the  characterization  of  the  Mie-Griineisen  EOS  for  isotropic  materials  by 
measuring  shock  velocity  Us  and  particle  velocity  up.  The  method  was  originally  employed  by 
Walsh  and  Christian  (91)  and  later  by  many  others  (e.g.,  McQueen  and  Marsh  [92];  Van  Thiel  et 
al.  [95];  Marsh  [94];  Steinberg  [49];  Meyers  [6];  Trunin  et  al.  [95]).  Using  Us-up  experimental 
data,  further  approximation  of  this  experimental  curve  using  equation  9  or  10  can  be  done.  The 
same  methodology  can  be  used  to  validate  the  modified  Mie-Griineisen  EOS  for  anisotropic 
materials;  thus,  a  numerical  simulation  of  the  anisotropic  plate  impact  test  is  considered. 

In  addition  to  its  fundamental  utility  in  characterizing  the  EOS,  the  plane  shock-wave  technique 
provides  a  powerful  tool  for  studying  material  properties  at  different  strain  rates  (e.g.,  Kipp  and 
Grady  [9(5];  Steinberg  [49];  Meyers  [6];  Johnson  et  al.  [50];  Kanel  et  al.  [57];  Millett  et  al.  [52]). 
Such  characteristics  as  spall  pressure,  shock  velocity,  particle  velocity,  HEL,  thickness  of  the  spall 
section,  time  to  spall,  and  free-surface  velocity  of  the  spall  section  can  be  measured  and  used  for 
the  characterization  of  material  dynamic  response  (e.g.,  Meyers  [(5];  Kanel  et  al.  [57];  Kanel  [45]; 
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Krivtsov  [97];  Millett  et  al.  [52];  Bourne  and  Gray  III  [6/  ];  and  Stoffel \9H\).  In  the  earliest 
plane- wave  experiments,  two  parameters  that  could  be  determined  were  the  HEL  (or  stress  level 
associated  with  the  elastic  precursor  wave)  and  the  dynamic  compressibility  (or  bulk  modulus) 
associated  with  the  following  plastic  wave.  More  complex  theories  of  material  behavior  have  led 
to  the  extraction  of  much  additional  information  from  such  experiments  (definition  of 
non-standard  parameters,  e.g.,  Kiselev  and  Lukyanov  [99]). 

From  an  experimental  point  of  view,  it  is  clear  that  the  production  of  waves  in  a  metal  target  and 
the  measurement  of  their  characteristics,  such  as  speed  and  intensity,  provide  one  of  the  most 
convenient  methods  of  investigating  the  physical  properties  of  a  material  under  high  pressure.  A 
schematic  of  the  plate  impact  test  is  depicted  in  figure  1,  in  which  a  flyer  plate  impacts  a  target 
plate,  which  is  bonded  to  a  polymethylmethacrylate  (PMMA)  plate.  In  this  report,  the  case  is 
considered  where  the  diameters  of  the  flyer  and  the  target  are  much  greater  than  their  thicknesses 
and  the  characteristic  time  of  the  process  is  the  time  of  several  runs  of  elastic  waves  across  the 
thickness  of  the  target  plate.  In  such  a  case,  the  problem  may  be  solved  using  a  uniaxial  strain 
state  (one-dimensional  mathematical  formulation  in  strain  space)  and  the  adiabatic 
approximation;  therefore,  planar  impact  generates  two  one-dimensional  shock  waves.  One 
propagates  into  the  target  and  the  other  into  the  flyer  plate.  These  shock  waves  reflect  as 
rarefaction  waves  from  the  free  surfaces  of  the  flyer  and  from  the  back  of  target  plates  connected 
to  the  PMMA.  With  a  thin  flyer,  these  rarefaction  waves  interact  inside  the  target,  producing  a 
state  of  tension  in  some  regions,  which  leads  to  the  spallation.  A  large  set  of  spall  experiments 
have  been  carried  out  by  Kanel  et  al.  (100,  101).  In  this  work,  the  dynamic  anisotropic 
elasto-plastic  material  behavior  before  spallation  is  considered. 

The  shock-wave  experiment  has  certain  potential  advantages  associated  with  the  level  of  strain 
rate  that  can  be  induced  and  the  commonly  accepted  belief  that  no  geometrical  dispersion  effects 
occur.  It  has  frequently  provided  the  motivation  for  the  construction  of  material  constitutive 
relations  and  has  been  the  principal  means  for  determining  material  parameters  for  some  of  these 
relations  (Steinberg  [49]\  Meyers  [6];  Bourne  and  Gray  III  [ 61]\  Chijioke  et  al.  \ 63 ] ;  Gu  and 
Ravichandran  [48]). 

1.2  Linear  Anisotropic  Constitutive  Relationship 

It  is  well  known  that  elastic  isotropic  materials  are  characterized  by  two  constants  (Lame 
parameters).  In  the  case  of  anisotropic  elastic  material,  the  stiffness  matrix  is  required  to  describe 
the  constitutive  relationship.  The  strain-stress  constitutive  relationship  for  any  elastic  material 
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Target  plate 


Gauge 


Figure  1.  Schematic  diagram  of  the  experimental  target  assembly. 


may  be  represented  in  contracted  form  as 


“U 


Jijkl&kl 


<7. 


Cijki^kl 


Cijkl 


-1 

ijkl 


(i) 


where  are  elements  of  the  stiffness  matrix  and  ,7ufc/  are  elements  of  the  compliance  matrix. 
In  this  report,  contraction  by  repeating  indexes  is  assumed.  In  general,  the  stiffness  matrix  Cijki 
may  be  a  function  of  aVJ ,  ,  etJ  and  some  other  parameters.  However,  it  is  somewhat  unwieldy 

as  such,  and  is  often  considered  to  be  constructed  of  constants,  which  produces  the  familiar 
Hooke’s  Law,  equation  1.  One  reason  why  the  deficiency  of  Hooke’s  Law  becomes  apparent 
experimentally  under  large  pressures  is  that  the  bulk  modulus  of  the  material  is  quite  different 
from  the  material’s  stress-free  value.  For  isotropic  materials,  this  problem  has  been  solved  by  the 
introduction  of  conventional  decomposition  into  two  quantities:  the  hydrostatic  stress  (or 
pressure),  which  only  induces  a  change  of  scale,  and  the  deviatoric  stress,  which  only  induces  a 
change  of  shape  (e.g.,  Wilkins  [102]).  This  detailed  description  of  such  well-known  facts  will  be 
used  later  in  the  construction  of  a  generalized  decomposition  of  the  stress  tensor,  which  will  take 
into  account  physical  properties  of  anisotropic  materials.  The  pressure  is  defined  as  one-third  the 
negative  sum  of  the  three  normal  stress  components.  For  isotropic  materials  in  the  elastic  region, 
the  pressure  is  directly  linked  with  volumetric  strain,  and  the  decomposition  of  the  stress  tensor 
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has  the  conventional  form 


Uij  =  -pSij  +  Sij  ,  p  =  —Ke  ,  Sij  =  2Gefj  ,  p  = -- (an  +  022  +  033)  ,  (2) 

where  p  is  the  hydrostatic  pressure,  SK]  is  the  deviatoric  part  of  the  stress  tensor,  K  is  the 
conventional  bulk  modulus,  G  is  the  shear  modulus,  d%]  is  the  Kronecker  delta  symbol, 
efj  =  6ij  —  \e8ij  is  the  deviatoric  strain  tensor,  and  e  =  en  +  e22  +  £33  is  the  volumetric  part  of  the 
strain  tensor.  The  distortional  change  is  represented  by  the  deviatoric  strain  tensor  efj.  Since 
experimental  evidence  reveals  that  the  compressibility  of  many  materials  changes  under  large 
pressures,  the  deviatoric  formulation  suggests  that  while  the  simplicity  of  Hooke’s  Law  (constant 
coefficients)  might  possibly  be  retained  for  computation  of  the  deviatoric  stresses  and  strains,  a 
more  accurate  scalar  EOS  should  simultaneously  be  employed  to  account  for  nonlinear 
compressibility  effects. 


2.  Isotropic  Shock- Wave  Modeling 


During  moderate  to  high  levels  of  shock  loading,  the  material  undergoes  nonlinear  behavior  in 
which  the  deformation  is  thermodynamically  coupled  with  the  internal  energy;  therefore,  an  EOS 
is  required  to  describe  the  material’s  response  to  these  conditions.  It  is  convenient  in  numerical 
codes  to  have  an  analytical  form  of  the  EOS,  but  such  an  analytic  form  is  at  best  an  approximation 
to  the  true  relationship.  The  EOS  for  the  computational  treatment  of  isotropic  materials  typically 
defines  the  pressure  as  a  function  of  density  p  (or  specific  volume,  v)  and  specific  internal  energy 
e.  A  very  popular  form  of  EOS  that  is  used  extensively  for  isotropic  solid  continua  is  the 
Mie-Griineisen  EOS: 

P  =  f{p,e)  =  pr{v)  +  (e  -  er(i/))  ,  (3) 

where  v  is  the  specific  volume  and  T(v)  is  the  Griineisen  gamma  defined  as 

.  (*)_  . 

Traditionally,  T  is  taken  to  be  constant  T  =  T0;  alternatively,  it  has  often  been  assumed  that 
Tq/vq  =  T(z/)/i/  =  constant  in  equation  3.  Functions  Pr{y )  and  er(y)  are  assumed  to  be  known 
functions  of  v  on  some  reference  curve.  Possible  reference  curves  include:  the  shock  Hugoniot 
curve,  a  standard  adiabatic  curve  (e.g.,  the  adiabatic  through  the  initial  state  (p0,  z/0)),  the  0  K 
isotherm,  the  isobar  p  —  0,  the  curve  e  =  0,  or  some  composite  curve  of  one  or  more  of  the  above 
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curves  to  cover  the  complete  range  of  interest  in  the  parameter  v.  The  most  commonly  used  form 
of  the  Mie-Griineisen  EOS  for  solid  materials,  which  uses  the  shock  Hugoniot  as  the  reference 
curve,  is  obtained  by  combining  equation  3  with  the  Rankine-Hugoniot  equations  to  obtain: 

P  =  f(p,e)  =  pH  ■  ^1  - +  pTe  ,  (5) 

where  PH  is  the  Hugoniot  pressure,  fi  =  p/p0  —  1  is  the  relative  change  of  volume,  T  is  the 
Griineisen  function,  p  is  the  density,  p0  is  the  initial  density,  e  is  the  specific  internal  energy,  and  v 
is  the  current  specific  volume.  The  shock  Hugoniot  is  often  preferred  for  the  Mie-Griineisen 
reference  curve,  precisely  because  the  experimental  data  used  to  fit  the  EOS  are  shock-Hugoniot 
data.  In  this  way,  the  Mie-Griineisen  reference  curve  need  not  be  derived  from  the  fitting  data, 
but  literally  represents  that  data  in  terms  of  the  Ph  function.  Furthermore,  there  is  a  perceived 
comfort  in  knowing  that  the  computational  trajectory  of  pressure  vs.  volume  for  an  impact 
simulations  will  lie  closer  to  the  known  reference  Hugoniot  curve  than,  for  example,  a  cold  curve. 

The  Rankine-Hugoniot  equations  for  the  shock  jump  conditions  can  be  regarded  as  defining  a 
relation  between  any  pair  of  the  p,  p,  e,  up  variables  and  Us  (Meyers \6\).  In  many  dynamic 
experiments,  up  (the  particle  velocity  directly  behind  the  shock)  and  U s  (the  velocity  at  which  the 
shock  wave  propagates  through  the  medium)  are  measured.  Generally,  the  shock  velocity  Us 
(and,  by  inference,  the  Hugoniot  pressure)  is  a  nonlinear  function  of  particle  velocity  up  and  it  is 
sometimes  fit  to  the  following  polynomial  relation  (Steinberg  [49))\ 

Us  —  C  +  SlUp  +  S2  (jjj  Up  +  S3  Up  .  (6) 

This  functional  form  is  more  often  employed  in  the  strictly  linear  form  (with  S2  and  S3  identically 
zero).  Kerley  (103)  shows  that  the  linear  Us-up  relationship,  because  of  its  wide  acceptance,  is 
often  treated  almost  as  a  “law”  rather  than  a  fit.  He  explains  that  the  general  linearity  observed  in 
the  Us-up  relationship  arises  because  it  is  basically  a  “plot  of  up  vs.  itself,”  and  that  the  linearity 
disappears  when  plotting  the  difference  Us  —  up  vs.  up .  Kerley ’s  notation  of  a  “plot  of  up 
vs.  itself”  refers  to  the  fact  that,  considering  Us  as  the  sum  of  up  and  (Us  —  up ),  the  magnitude  of 
(Us  —  up)  is  significantly  smaller  than  that  of  up  for  strong  shock.  Since  up  vs.  up  is 
tautologically  linear,  the  term  (Us  —  up)  therefore  carries  more  concentrated  EOS  information 
than  Us  alone,  assuming  one  can  measure  it  to  as  comparable  an  accuracy  as  Us.  Nonetheless, 
most  shock-Hugoniot  data  fits  are  provided  in  this  linear  Us  vs.  up  form. 

According  to  the  constraints  of  Griineisen  theory,  the  Griineisen  parameter  T  must  be  a  function 
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of  volume  alone.  Many  empirical  forms  have  been  proposed  including 


7o  +  ap 
1  +  /I 


(7) 


This  form  reduces  to  the  simpler  fit,  namely,  u/T  =  constant,  mentioned  earlier,  for  the  case 
where  a  is  zero.  One  must,  however,  be  mindful  of  the  compatibility  of  fits  between  the 
Hugoniot  and  Griineisen  functions.  Segletes  ( 104-106 )  has  shown  that  arbitrary  specification  of 
shock-Hugoniot  fits  and  Griineisen  fits  can  lead  to  thermodynamic  instability  of  the  EOS  model. 


Nonetheless,  the  Griineisen  EOS,  if  one  adopts  the  cubic  shock-velocity  relation  of  equation  6, 
can  be  derived  as 


P  =  T 
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(8) 


when  p  >  0,  and 


p  =  Po<?p  +  (1  +  p)  ■  T  •  E 


(9) 


when  p  <  0.  In  these  equations,  E  is  the  internal  energy  per  initial  specific  volume,  c  is  the 
intercept  of  the  Us~up  curve  (nominally  akin  to  the  ambient  bulk  sound  speed  of  the  material,  c0), 
Si,  S2,  and  S:i  are  the  coefficients  of  the  slope  of  the  Us-up  curve  equation  6,  y0  is  the  ambient 
Griineisen  gamma,  and  a  is  the  first  order  volume  correction  to  70.  Parameters  c,  Si,  S2,  S3,  70, 
and  a  represent  material  properties,  which  define  its  EOS.  Parameters  have  been  defined  to  cover 
a  large  number  of  isotropic  materials  (Steinberg  [49]).  The  tensile  part  of  the  EOS,  equation  9,  is 
outside  the  domain  of  the  shock  Hugoniot,  but  instead  assumes  a  constant  bulk  modulus  in  the 
hydrostatic  portion  of  the  equation,  along  with  Griineisen  energy  coupling. 


Other  forms  of  the  Mie-Griineisen  EOS  have  been  put  forth,  based  upon  different  assumptions 
regarding  the  reference  curve  and  Griineisen  function,  for  example,  a  frequency-based  approach 
(Segletes  and  Walters  [107]\  Segletes  [108,  109]),  in  which  the  Mie-Griineisen  reference  curve 
and  the  Griineisen  function  are  not  independently  specified,  but  rather  are  both  interdependent 
functions  of  the  characteristic  frequency  of  the  lattice,  u>: 

Cp  \ 

r  0k) 

where  the  subscript  0  refers  to  a  parameter  at  the  ambient  state,  k  is  a  constant  parameter  defined 


{  [{u/uoY  -  1]  +  k  (k  -  1)  (uj/ujoT  In  (cu/co’o)}  ,  (10) 


pu/T  —  E  = 
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by 


ft  —  [rvoio  —  1  +  3/2  ■  ro£]/To  , 


(11) 


c  is  the  bulk  sound  speed,  Tvo]  is  the  so-called  “volumetric”  Griineisen  parameter  (/.<?.,  the 
Dugdale-MacDonald  [110]  definition),  and  the  constant  £,  while  defined  by 


£  =  4/(3r0)  -  d(u/T)/du\ 0 


(12) 


is  generally  fit  directly  to  EOS  data  (since  dT/du  is  not  easily  measured). 


Griineisen  theory  indicates  that  the  characteristic  lattice  vibrational  frequency  is  related  to  the 
Griineisen  parameter  by  way  of  T/v  =  ( dp/dE)u  =  — ( duj/dv)/u .  Thus,  by  employing  a 
preferred  functional  form  for  the  volume-dependent  lattice  frequency, 


cc/cco 


(z/0/z/)i/3  -  1 


(13) 


the  volume-dependent  function  for  T  may  be  obtained  directly  through  differentiation  of  u, 
thereby  completely  defining  the  EOS  of  equation  10  operationally  in  terms  of  p,  E,  and  u,  given 
the  five  material  constants  defining  the  columns  of  table  1. 


Table  1.  EOS  parameters  for  select  materials. 


C0  (m/s) 

l/v0  (kg/m3) 

r0 

EvolO 

£ 

Ag 

3221 

10490 

2.22 

2.29 

0.500 

A1 

5189 

2700 

2.03 

1.84 

0.716 

Cu 

3995 

8930 

2.02 

2.10 

0.520 

Stainless  Steel 

4571 

7896 

1.81 

2.00 

0.550 

From  the  general  EOS  form  of  equation  10,  expressions  may  be  derived  for  the  cold  curve,  pc, 

Pc=  (Co /To)2  (cj/uj0)k  In  (uj/ujo)/^  ,  (14) 


as  well  as  the  Hugoniot, 

CqV 

t0k/ 

In  both  of  these  equations,  the  term  t/)  =  v/T  is  used  as  convenient  shorthand.  These  curves  may 
be  compared  to  both  diamond-anvil  (essentially  “cold”  compression)  and  shock-Hugoniot  data. 


{  [(u/ujoT  -  1]  +  «(«  -  1)  (cc/o;o)Kln(a;/a;o)}  .  (15) 


Ph  [i>  ~  (vo  -  v)/2]  = 
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An  example  of  such  a  comparison  to  aluminum,  for  two  different  levels  of  zoom,  is  given  in 
figure  2.  As  can  be  seen  from  the  nature  of  equations  14  and  15,  there  is  no  intrinsic  way, 
through  the  selection  of  the  parameters,  to  independently  fit  both  cold  and  Hugoniot  curves.  That 
both  curves,  nonetheless,  fit  the  data  so  well  (and  for  a  variety  of  different  materials)  is  indicative 
of  a  more  general  suitability  of  the  parent  form  described  by  equation  10. 


Figure  2.  Cold-compression  and  shock-Hugoniot  curves  for  aluminum  (a)  to  2.4  and  (b)  1 1  megabars. 
Note  that  cold-compression  data  arc  filled  symbols  and  Flugoniot  data  arc  open  symbols. 
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3.  Anisotropic  Shock- Wave  Modeling 


Various  models  and  numerical  methods  have  been  developed  to  simulate  processes  in  solid  media. 
The  choice  of  these  methods  will  depend  largely  on  the  problem  and  the  available  resources.  A 
detailed  review  of  the  most  common  of  these  was  provided  by  Benson  (111). 

3.1  Characteristics-based  Solver 

For  modeling  physical  problems,  such  as  high  velocity  impact  and  explosive  loading  of  isotropic 
and  anisotropic  solid  materials,  there  is  strong  demand  on  the  numerical  methods  to 
simultaneously  achieve  high  shock-wave  resolution,  maintain  sharp  interfaces,  and  accurately 
impose  interfacial  boundary  conditions  subject  to  finite  deformation.  A  characteristics-based 
solver  can  be  successfully  used  for  solving  the  resulting  hyperbolic  partial  differential  equations 
(PDEs).  Depending  on  the  constitutive  assumptions  regarding  the  stress-strain  relation,  the 
resulting  hyperbolic  PDEs  will  have  a  different  form  (i. e. ,  different  Lagrangian).  Conservative 
formulations  of  the  governing  laws  (PDEs)  of  elasto-plastic  solid  media  have  distinct  advantages 
when  solved  using  high-order  shock  capturing  methods  for  simulating  processes  involving  large 
deformations  and  shock  waves.  Let  us  define  the  density  of  Lagrangian  L  for  our  media,  per  unit 
initial  volume,  as 

L  =  ift,  (X)  X  x\  -  E  (x,  S) 


E  (x,  S)  —  E  (X,  iifj,  S) 
E  (x,  S)  =  p0e  (x,  S) 


(16) 


^  =  \  (£  FuFekj  -  % 

where  rf-}  =  rfi3  (X,  t)  are  the  elastic  Lagrangian  strains  from  X  to  x,  r/?;j  =  rpj  (X,  t)  are  the  total 
Lagrangian  strains  from  X  to  x,  LjC  are  elements  of  the  elastic  deformation  gradient,  Fkj  are 
elements  of  the  deformation  gradient,  and  E  (x,  S)  is  the  density  of  internal  energy  per  unit  initial 
volume.  Lagrange’s  equations  of  motion  in  the  absence  of  body  forces  are 


d  d  L 
dt  dii 


d  dL 

dXkWk 


0 


(17) 
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or 


(18) 


Note  that 

^  l  (FaStm  +  FimS„i)  ,  =  +  ■  (W) 

o^ik  Z  Orik  Z 

With  equation  19,  and  noting  the  symmetry  of  r)t]  and  rf^,  the  equation  of  motion  equation  18  is 


Po  (X)ii 


(20) 


A  characteristics-based  solver  can  be  applied  when  the  system  equation  20  is  closed  by  analytical 
formulae  for  the  specific  internal  energy  e  (x,  S).  For  small-amplitude  elastic  anisotropic  wave 
propagation  (rjfj  =  r)l3),  specific  internal  energy  can  be  written  as  quadratic  form: 

pe  (x,  S )  =  +  TS  .  (21) 

For  isothermal  motion,  the  potential  in  the  Lagrangian  density  L  is  the  density  of  free  energy 
F  (x,  S)  =  F  (X,  rfi } ,  S)  per  unit  initial  density  in  place  of  E  (x,  S)  =  E  (X,  rjfj,  S ) ,  and  the 
adiabatic  elastic  constants  C^kl  are  then  isothermal  Cfjkl.  In  the  case  of  isotropy,  an  isentropic 
hyperelastic  EOS  in  terms  of  the  invariants  of  the  elastic  Greens  tensor  was  considered  (Miller 
and  Colella  [112]:  Barton  et  al.  [113]): 

e  M,  j3)  =  ~  p  M  +  F  (ji  -  3 4/3) 


Ji  —  tr  (C) 

J2  =  \  [(trC)2  -trC2] 


(22) 


J3  =  det  |C| 

where  C  =  (F'  )7  Fe  is  the  elastic  Green’s  tensor.  The  quantity  P  (y)  can  be  described  by  any 
form  of  isotropic  EOS.  Another  form  of  the  isotropic  hyperelastic  EOS  has  been  put  forth 
(Dorovskii  et  al.  [114]:  Barton  et  al.  [113]:  Barton  and  Drikakis  [11 5]): 

e  (Ju  J2,  J3,  S)  =  ^  (j3“/2  -  1  y  +  CuT0Jj/2  (exp  [S/cv]  -  l)  +  =yjf 2  (J2/ 3  -  J2)  ,  (23) 
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where  J\  =  tr  (Cf),  J2  =  \  [(tr  CF)2  —  tr  C|>] ,  and  J3  =  det  \CF\  are  the  invariants 
corresponding  to  the  elastic  Finger  tensor  CF  —  (Fe)~r(Fe)-1.  The  parameters 
K0  =  Cq  —  (|)  Bo  =  bo  are  the  squared  bulk  speed  of  sound  and  the  squared  speed  of  shear 
waves,  respectively;  cv  is  the  heat  capacity  at  constant  volume;  and  a,  3,  7  are  constants 
characterizing  the  nonlinear  dependence  of  sound  speeds  and  temperature  on  the  mass  density. 
The  specific  internal  energy  can  be  decomposed  into  potentials  describing  the  cold  compression, 
ec  ip,  S);  thermal  energy  density,  e,  (p,  S);  and  the  contribution  due  to  shear  strain,  es  (Fe,  S).  A 
more  general  function  might  also  include  a  contribution  due  to  changes  in  dislocation  density, 
eh  (Fe,  S)  (Schreyer  and  Maudlin  [116]).  Thus,  a  general  form  could  be  written 

e  (Fe,  S)  =  ec  (p,  S )  +  e,  (p,  S)  +  e,  (Fe,  5)  +  eh  (Fe,  S)  .  (24) 


High-order  shock  capturing  methods,  based  upon  solving  Riemann  problems  locally  at  each  cell 
edge  throughout  a  computational  domain,  have  emerged  as  a  favorable  approach  to  meeting  the 
shock-wave  modeling  requirements.  Recently,  a  high-order  shock  capturing  scheme  was 
proposed  for  solid  dynamics  where  Godunov’s  method  was  applied  to  an  Eulerian  model 
(Godunov  and  Romcnskii  [/ /7J)  in  conservative  form  using  fixed  Cartesian  grids  (Miller  and 
Colella  [ 112];  Barton  et  al.  [113];  Barton  and  Drikakis  [ii5]).  Application  of  these  numerical 
methods  to  solid  mechanics  has  been  made  possible  by  formulations  of  the  governing  theory  as 
first-order  hyperbolic  systems  of  conservation  laws  in  the  Eulerian  frame  (Godunov  and 
Romenskii  [117];  Kondaurov  [118];  Plohr  and  Sharp  [119];  Godunov  and  Romenski  [120]). 
Introducing  the  vector  of  primitive  variables  W,  equation  20  can  be  rewritten  as  a  quasi-linear 
system  naturally  suitable  for  a  characteristics-based  solver: 


dW 

~dt 


A' 


tdW 

dxa 


Qp 


(25) 


where  Qp  is  the  vector  related  to  irreversible  processes  within  the  system.  In  the  ensuing 
computational  method,  the  convective  flux  terms  in  equation  20  are  usually  discretized  using  the 
well-known  method  of  Godunov.  The  solution  is  therefore  required  of  a  local  Riemann  problem 
at  the  boundaries  of  each  cell  in  the  computational  mesh.  The  solution  is  found  using  an 
approximate  method  based  upon  the  characteristic  tracing  and  thus  requires  detailed  knowledge 
of  the  eigen-values  and  the  eigenvectors  of  equation  20. 


The  method  of  characteristics  was  successfully  used  to  demonstrate  the  accuracy  of  the 
lamination  theory  (Whitney  and  Pagano  [121])  in  transient  wave  propagation  problems  for 
anisotropic  materials.  This  method  is  adequate  for  stiff  materials  or  weak  shocks  where  the 
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compression  is  small  and  the  elastic  response  is  nearly  linear.  The  laminated  plate  equations  used 
in  this  report  are  those  derived  by  Whitney  and  Pagano  (121),  which  are  considered  as  a 
representative  “first-order”  laminated  plate  theory.  Note  that  several  analytical  and  numerical 
solutions  to  the  laminated  plate  equations  were  constructed  in  the  past.  For  example,  Wang  and 
Tuckmantal  (722),  and  Moon  (/ 2.3)  have  discussed  the  wave  surface  based  on  the  laminated  plate 
equations  due  to  abrupt  change  in  stress.  Also,  Chow  (124)  applied  the  Laplace  transform 
technique  and  calculated  the  deflection  of  a  laminated  plate  under  a  concentrated  load. 

Moon  (123)  solved  the  problem  of  one-dimensional  wave  propagation  due  to  a  prescribed  line 
load  by  the  fast  Fourier  transform  technique.  Furthermore,  Sun  and  Lai  (125)  analyzed  the 
response  of  a  unidirectional  fiber-reinforced  layer  subjected  the  response  of  a  unidirectional 
fiber-reinforced  layer  subject  to  the  laminated  plate  equations  and  exact  equations.  For 
symmetric  cross-ply,  or  balanced  angle-ply,  laminated  plates  under  one-dimensional  in  plate 
impact,  the  resulting  hyperbolic  partial  differential  equations  (Whitney  and  Pagano  [727])  reduce 
to  (Mortimer  et  al.  [69]) 

.  (26) 

where  x  represents  the  spatial  coordinate  in  the  direction  of  wave  propagation,  t  the  time 
coordinate,  u°  (x,  t)  the  mid-plane  displacements  in  the  x-direction,  and  M the  plate  mass  term. 
This  is  a  simple  wave  equation  with  wave  velocity  c0  =  (/1|  \/M) l/2.  For  a  finite  length  striker 
plate  of  the  same  material  as  the  specimen,  the  response  in  the  specimen  is  simply  a  rectangular 
wave  of  amplitude  e0  =  V0/(2c0)  and  pulse  length  t0  =  2//c0,  where  e0  is  the  strain,  V0  the  initial 
velocity  of  the  striker,  and  l  is  the  length  of  the  striker. 

For  symmetric  cross-ply,  or  balance  angle-ply,  laminated  plates  under  one-dimensional  transverse 
shear-bending  impact,  the  governing  equations  (without  surface  tractions)  reduce  to 

(Ain 

kA55  (Wl  +  v2tu)  =  M—  ,  (27) 

Dn'V2‘ipx  +  -DieV2^  —  kAss  (ijjx  +  Vic)  =  /  ,  (28) 

DieV2ripx  +  Dqq V2^  —  kA^i/jy  =  I  ,  (29) 

where  the  Al3,  Di;),  and  /  are  plate  coefficients  as  defined  in  Whitney  and  Pagano  (121);  w,  V'x,  U'y 
the  specimen  deflection  and  rotations,  respectively;  and  k  the  shear  correction  factor.  These 
equations  27-29  are  a  set  of  totally  hyperbolic  partial  differential  equations  and  were  successfully 
solved  by  a  numerical  method  of  characteristics  by  Mortimer  et  al.  (69). 
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3.2  Nonlinear  Anisotropic  Elastic  Incremental  Formalism 


In  the  early  1970s,  Johnson  et  al.  (31-33)  developed  a  continuum  framework  for  describing  the 
rate  dependent  elastic-plastic  response  of  single  crystals  to  shock-wave  loading.  Their  tensor 
formulation  coupled  linear  elastic  response  to  a  dislocation  dynamics  model  that  incorporated  the 
appropriate  slip  system  in  the  shocked  crystals.  In  this  framework,  both  elastic  and  plastic 
anisotropy  were  explicitly  accounted  for.  This  approach,  similar  to  the  above  method  of 
characteristics,  is  adequate  for  stiff  materials  or  weak  shocks,  where  the  compression  is  small  and 
the  elastic  response  is  nearly  linear.  However,  for  more  compliant  materials,  crystals,  or  stronger 
shocks,  several  factors  must  be  accounted  for:  (1)  large  compressions  lead  to  a  nonlinear  elastic 
response,  (2)  low  yield  stresses  lead  to  large  deformations,  and  (3)  large  compressions  and  large 
plastic  deformations  lead  to  a  significant  temperature  increase.  As  a  result,  to  address  these 
complications,  Winey  and  Gupta  (28-30)  developed  a  consistent  Lagrangian  thermomechanical 
framework  for  modeling  the  response  of  a  single  crystal  under  shock  loading.  This  framework  is 
also  applicable  for  anisotropic  materials,  highly  anisotropic  materials  such  as  molecular  crystals, 
and  strong  shocks  in  other  single  crystals.  The  Winey  and  Gupta (28-30)  approach  is  constructed 
on  the  thermodynamic  approach  of  Wallace  (126,  127),  which  couples  nonlinear  elasticity  within 
a  thermodynamically  consistent  incremental  formalism.  In  this  method,  the  incremental 
displacement  from  some  intermediate  configuration  xf  at  time  n  to  the  next  configuration  a;"'f  1  at 
time  n+1, 


(Ui 


n+1  —  r«+l  _ 

n  i 


(30) 


where  (tq)™+1  denotes  the  displacement  at  time  n  +  1  relative  to  time  n.  The  incremental  strains 
(negative  in  compression),  given  by 


Aey 


1 

2 


flfror1 1  9mt1 

dx1-  dx" 


(31) 


are  measured  relative  to  the  configuration  at  time  n.  The  reference  configuration  is  updated  from 
the  current  configuration  to  the  next  configuration  after  every  strain  increment.  The  incremental 
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equations  for  internal  energy  (Gibbs  equality)  are 


En+ 1  —  En  =  —  +7)n+1/2  Aei7-  +  Tn+1/2AS 

pTl  '  J  '  71  J 

(*«):+1/2  =  \  i+r1 + (%>:]  .  <32> 


j~in+ 1/2  _  _  ^n+l  _|_  rjin' 


where  pn  is  the  density  for  the  reference  configuration  at  time  n;  Tn+1  and  Tn  are  the  temperature 
of  the  configuration  at  time  n+1  and  n,  respectively;  and  ++'  ’  1  and  (+■)"  are  the 
thermodynamic  stresses  for  configuration  at  time  n  +  1  and  n,  respectively.  The  thermodynamic 
definitions  of  temperature  and  stress  tensor  are  defined  by 


rpn 


dEn 

~dS 


eij 


(tij  ) 


n+1 

n 


p 


d  En+1 


(33) 


Similar  to  Wallace  (126),  the  anisotropic  Griineisen  tensor  was  defined  by  Winey  and 
Gupta  (28-30): 


(r+:  = 


pnTr 


d(t 


n+1 
l3  )  n 


dS 


(34) 


J  s 


Therefore,  the  incremental  change  in  the  thermodynamic  stresses,  using  the  EOS  equation  34,  is 


(*«);+1  -  (««>:  =  (c«+r1/2  A  £«  -  r  (r„);+1/2  T^AS 


(35) 


where  the  linear  elastic  coefficients  are  defined  by 


(C, 


ijkl) 


d(t. 


lJ  /  n 


9e 


fcZ 


(36) 


In  the  Winey  and  Gupta  (28)  approach,  the  elastic  coefficients  in  equation  36  are  functions  of 
strain  and  entropy.  As  a  result,  the  incremental  change  from  configuration  at  time  n  to 
configuration  at  time  n  +  1  is 


(cijki)  r1  -  (cm)nn = 


oa. 


ijkl 


deri 


1  n+l/2 


Adm.n  T 


sa 


ijkl 


dS 


n+l/2 


•  AS 


(37) 


(Cijklmn)nn  = 


dC,, 


de 


ijkl 

mn  / 
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where  nonlinear  elastic  coefficients  CV,^mn  are  assumed  to  be  independent  of  strain  and  entropy 
in  this  formulation,  so  that  (C'lJ^mn)n+1  =  ( Cijkimn)n ■  To  maintain  thermodynamic  consistency, 
the  Winey  and  Gupta  (28)  approach  assumes  for  derivatives  that 


dCijki 

dS 


eij 


~P 


3(TTy) 


dej-i 


s 


(38) 


In  general,  is  a  function  of  entropy  and  strain.  However,  for  simplicity,  the  Winey  and 
Gupta  (28)  approach  assumes  T;3  to  be  constant.  Hence, 


dCijki 

dS 


eij 


-pT 


*? 


dT 


de 


ki 


=  pTYijTki 


(39) 


When  equation  39  is  substituted  into  equation  37,  the  incremental  change  in  the  linear  elastic 
coefficients  under  the  Winey  and  Gupta  (28)  assumptions  can  be  written  as 


( n  '\n+1  _  ( n  \n  —  (n  \n+l/2  A  ,  nrpn+ 1/2  /p  \n+ 1/2  /-p  \n+ 1/2  a  q 

y^ijkljn  y^ijkl) n  y^ijklmn) n  *  mn  '  P  ij)n  \  kl)n  * 


(40) 


To  complete  the  thermodynamic  description,  the  temperature  increment  from  configuration  at 
time  n  to  configuration  at  time  n  +  1  is  given  by 

(rri  \  77.+  1/2 

-J  AS  ,  (41) 

where  ce  is  the  specific  heat  capacity  at  constant  elastic  configuration.  To  maintain  consistency 
with  constant  Vl3 ,  it  is  assumed  that  ce  is  to  be  constant  as  well.  It  is  important  to  note  that  in  this 
formulation,  the  thermodynamic  stresses,  elastic  coefficients,  and  other  tensor  properties  must  be 
updated  during  the  incremental  calculation  so  that  they  are  referred  to  the  configuration  at  time  n. 

Using  the  above  framework,  the  propagation  of  large  amplitude  stress  waves  along  the  arbitrary 
directions  in  quartz,  sapphire,  LiF,  copper  single  crystal,  and  unreacted  PETN  are  performed  by 
Winey  and  Gupta  (28-30).  Measured  (Asay  el  al.  [73])  and  simulated  (Winey  and  Gupta  [29]) 
longitudinal  stress  histories  for  LiF  single  crystals  shocked  along  the  (100)  orientation  are 
presented  in  figure  3.  They  have  also  discussed  differences  between  pure  mode  wave  propagation 
for  linear  and  nonlinear  elastic  deformation. 


3.3  Conventional  Decomposition  of  Anisotropic  Equations 

Speaking  generally,  any  second-order  tensor  can  be  decomposed  into  the  spherical  part  and  the 
deviatorical  part.  In  the  case  of  continuum  mechanics,  the  decomposition  of  the  stress  tensor  and 
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Stress  at  Gauge  Interface  (GPa) 


Figure  3.  Measured  (solid  lines)  and  simulated  (dashed  lines)  longitudinal  stress  histories  for 
LiF  single  crystals  shocked  along  the  (100)  orientation.  The  experimental  data 
(Asay  et  al.)(73)  and  simulations  (Winey  and  Gupta) (29)  arc  for  a  6061-T6 
aluminum  plate  impacting  the  LiF  crystals  at  a  velocity  of  ca.  0.34  km/s,  where  the 
target  crystals  are  backed  by  a  quartz  stress  gauge.  The  stress  is  measured  at  the 
interface  between  the  LiF  crystal  and  the  quartz  gauge.  The  numbers  above  the 
curves  indicate  the  sample  thickness  in  mm.  Time  is  relative  to  the  moment  of 
impact.  Reprinted  with  permission  from  Winey,  J.M.,  Gupta,  Y.  M., 

J.  Appl.  Phys.  99,  023510,  (2006).  Copyright  2006,  American  Institute  of  Physics. 


19 


the  strain  tensor  into  their  volumetric  and  deviatoric  components  ( e.g .,  Wilkins  [702])  has  certain 
physical  justification.  This  step  is  done  in  order  to  distinguish  between  thermodynamic  (EOS) 
response  and  the  ability  of  the  material  to  carry  shear  loads  (strength).  For  anisotropic  materials, 
the  decomposition  of  the  stress  and  strain  tensors  into  spherical  and  deviatoric  parts  in  stress 
space  and  strain  space  results  in  stress  and  strain  components  that  do  not  correspond  to  each  other 
due  to  the  material  properties’  anisotropy.  For  example,  using  the  stress-strain  relation  for  an 
orthotropic  material,  the  hydrostatic  pressure  can  be  rewritten  as 


p  =  —  {£)(  1  -  u23u32)  +  E2(l  -  Z/13Z/31)  +  E3(  1  -  v12v2i) 

+2  ■  [E\  (/y2i  +  22 31^23)  +  Ei(u3  1  +  v2iv32)  +  E2(v32  +  ^31^12)]}  •  e 

—  ^  —  ^23^32)  +  Ei(v2i  +  v3iv23)  +  Ei(u31  +  2221z/32)}  • 

—  7^  {E\  {y-2\  +  223iZ/23)  +  E2(  1  —  2213z/31)  +  E2(v32  +  z/12z/31)}  •  e22 

—  {El  (2231  +  22212232)  +  E2  ( Z232  +  22i22231)  +  E3(l  —  2212Z221)}  •  633 


where  /3  =  1  -  2212z221  -  u13u31  -  v23v32  -  2u21u32u13,  Eu  E2,  E3  are  the  Young’s  moduli,  and 
are  the  six  Poisson  ratios.  Hence,  for  an  anisotropic  material  in  general,  the  mean  stress  depends 
on  the  deviatoric  strains  and,  as  a  result,  the  decomposition  used  for  isotropic  materials  is  not 
applicable. 


3.4  Deviatoric  Constitutive  Relationship  for  Anisotropic  Materials 


While  the  mathematics  and  physics  of  the  constant  coefficient  constitutive  relationship  for 
isotropic  materials  was  well  understood,  the  casting  of  these  rules  into  an  anisotropic  framework 
was  not  a  straightforward  task.  In  particular,  the  difficulties  were  associated  with  two  primary 
differences  in  the  behavior  of  anisotropic  materials  with  respect  to  that  of  isotropic  materials: 

(1)  under  hydrostatic  pressure,  strain  is  not  uniform  in  all  three  directions  of  the  material 
coordinates,  and  (2)  except  under  restrictive  modulus  conditions,  deviatoric  stress  will  produce 
volumetric  dilatation.  As  a  solution  of  these  difficulties,  the  elastic  generalized  deviatoric 
anisotropy  decomposition  was  proposed  by  Segletes  (34).  While  his  arguments  were  applied  to 
the  case  of  transversely  isotropic  materials,  the  generalization  to  orthotropy  is  straightforward  and 
is  shown  directly.  Decomposition  of  the  stress  and  strain  tensors  into  their  hydrostatic  and 
generalized  deviatoric  components  yields  (Segletes  [34)) 

Sij  =  < Tij  —  &5ij  ,  SijSij  =  0  ,  a  =  —  (an  +  <t22  +  cr33)  ,  (43) 

c  ij  2  tj  f.jj  ,  (44) 
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(45) 
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=  e  7^  0 


(46) 


-<?.j  ^ ijkl  -'Id 

where  a  is  the  hydrostatic  stress,  represents  the  normal  strains  due  to  hydrostatic  stress,  e  is  the 
dilation  of  generalized  deviatoric  strain,  and  e  is  the  excess  anisotropic  dilation.  One  may  acquire 
upon  substitution  into  equations  43-45 


Sij  4"  crSij  Cijki^ki  +  CijkiEki  .  (47) 

Unlike  the  isotropic  materials  in  which  a  hydrostatic  pressure  produces  a  uniform  dilatation  in  all 
three  coordinate  directions,  hydrostatic  strain  for  an  anisotropic  material  is  non-uniform. 
Therefore,  if  one  defines  the  deviatoric  components  of  stress  and  strain  to  be  the  total  stress  / 
strain  components  decremented  by  an  amount  that  would  result  from  a  hydrostatic  stress  state, 
one  can  conclude  (per  condition  (1)  above)  that  there  is  a  unique  hydrostatic  strain  component 
associated  with  all  three  directions  in  the  material  coordinates  (the  coordinate  system  that 
produces  no  shear  coupling).  Equation  47  may  be  decoupled  to  give  a  hydrostatic  equation, 


ij  Cijkl^kl  i 


(48) 


and  a  deviatoric  relationship  void  of  hydrostatic  terms 


Sij  Cijkieki 


(49) 


Under  the  influence  of  a  purely  hydrostatic  stress  state  (and  assuming  the  moduli  to  be  constant), 
there  will  be  constant  ratios  between  the  components  of  normal  strains  etj  due  to  hydrostatic 
stress.  Defining  the  ratios  in  terms  of  material  compliances  Jijku  it  follows  from  equation  45  that 


Cll 

^22 


K{2  = 


•hill  +  Jll22  +  J 1133 
J2211  +  J2222  +  <72233 


^33 

C22 


7^2  = 


<73311  +  J3322  +  <73333 
<72211  +  J2222  +  <72233 


(50) 


Recall  from  equation  46  that  the  sum  of  the  three  normal  deviatoric  strain  increments  is  not 
generally  zero,  but  rather  equals  a  deviatoric  dilatation,  e  .  The  significance  of  this  term  is  that  a 
state  of  stress  whose  average  normal  value  is  zero  can  produce  volumetric  change  on  an  element 
with  respect  to  that  element’s  stress  free  volume.  We  may  substitute  the  normal  components  of 
equation  44  into  equation  46  and,  with  the  aid  of  equation  50,  develop  a  relation  between  e  and 


^22 ' 


(en  +  e22  +  ^33)  —  e 
7^12  +  1  +  7t<32 


(51) 
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Likewise,  one  may  employ  the  condition  of  uniform  strain  (en  =  e22  =  e33)  into  the  deviatoric 
constitutive  relation,  equation  49,  to  acquire  constant  ratios  between  the  components  of  normal 
stresses  resulting  from  the  uniform  strain  state: 

°~11  _  _  Cllll  +  C 1122  +  Cll33  0_ 33  _  _  ^3311  +  C3322  +  C*3333 

&22  12  CWl  +  C*2222  +  C*2233  °22  ^  C'22 1 1  +  6-2222  +  C*2233 

So  as  to  insure  that  the  deviatoric  stress  has  no  hydrostatic  component,  substitute  equation  49  into 
equation  43,  in  light  of  equation  52,  to  show  that 


-^12ell  +  e22  +  -^32e33  —  0  .  (53) 

From  this  equation,  en  +  e22  +  e33  may  be  isolated  and  replaced  with  e,  and  the  remaining 
terms  eliminated  via  equation  44,  eventually  producing  a  second  relation  between  e  and  : 

-  6  +  (-ft  12  ~~  l)eH  +  (-ft 32  ~~  l)e33  <-4, 

22  K‘12(K?2  - 1)  +  iq2(iq2  - 1)  ' 

The  e22  term  may  be  eliminated  between  equations  5 1  and  54  to  produce  a  closed-form  for  e  in 
terms  of  the  current  strain  state  f  and  the  material  parameters  K"  and  K^.  Once  e  is  obtained, 
e22  follows  from  either  equation  5 1  or  54.  The  other  f  r]  components  follow  from  equation  50, 
and  finally,  from  equation  44,  the  deviatoric  strain  terms  follow.  All  unknown  deviatoric 
strains  are  now  solved.  As  an  aside,  to  recreate  the  transversely  isotropic  condition  described  by 
Segletes  (34),  the  terms  I\Z2  and  KZ2  need  merely  be  set  equal  to  unity. 

Summing  the  three  equations  for  normal  stress  equation  48  subject  to  equations  44  and  46  yields 
upon  reduction 

cr  =  Ka  (e 1 1  +  622  +  £33  +  e)  ,  (55) 

where  Ka  is  a  true  material  property,  which  is  called  the  first  effective  bulk  modulus  of  the 
material  (Lomakin  [128]).  This  modulus  equals  one  ninth  the  sum  of  the  nine  normal  stiffness 
matrix  components  67^: 

1  1  3  3 

Ka  =  -SjjCjjkiSki  —  -  Cjj  ,  (56) 

2=1  j 

where  Cl3  are  elements  of  the  stiffness  matrix  (written  in  Voigt  notation). 

Alternately,  the  use  of  the  deviatoric  constitutive  relation,  equation  49,  hinged  upon  the 
satisfaction  of  equation  48.  Inverting  equation  49  into  compliance  form  and  summing  the  three 
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equations  for  normal  strain  yields  upon  reduction 


a  =  Ke  (en  +  622  +  e33  -  e) 


(57) 


where  Ke  is  a  true  material  property,  which  is  called  the  second  effective  bulk  modulus  of  the 
material  (Segletes  [34]).  This  modulus  equals  the  reciprocal  of  the  sum  of  the  nine  normal 
compliance  matrix  components  Jijki'. 


6{j  'J/jkid/j 


1 

~ 3  3 

E  E  -hj 

i=l  j= 1 


(58) 


where  Jt]  are  elements  of  the  compliance  matrix  (written  in  Voigt  notation).  These  effective 
moduli,  unlike  the  bulk  modulus,  are  independent  of  deviatoric  stress  in  anisotropic  materials. 

The  first  and  second  bulk  moduli  reduce  to  the  conventional  bulk  modulus  in  the  limit  of  isotropy. 
Once  a  is  known,  the  deviatoric  stresses  follow  from  equation  43. 


It  was  mentioned  previously  that,  for  real  materials  under  large  compression,  the  empirical 
relation  between  dilatation  and  pressure  is  not  a  linear  one  and  at  high  pressures  is 
thermodynamically  coupled  to  internal  energy.  One  advantage  of  this  deviatoric  formulation  lies 
in  the  ability  to  arbitrarily  make  the  hydrostatic  relations  equations  57  and  58  nonlinear  while 
retaining  the  linear  simplicity  of  Hooke’s  Law  for  the  deviatoric  portion  of  the  constitutive 
relation.  Though  this  ad  hoc  procedure  does  not  theoretically  follow  as  an  extension  to  Hooke’s 
Law,  it  does  permit  the  code  user  to  more  flexibly  model  the  empirical  EOS  behavior  of  the 
material,  as  —  a  =  p  (y/ exp(e),  E).  These  tensors  differ  from  the  absolute  stress  /  strain  tensors 
in  that  the  normal  components  of  stress  and  strain  are  decremented  by  the  hydrostatic  values  of 
the  normal  stresses  and  strains,  respectively.  In  this  way,  the  deviatoric  quantities  represent 
deviation  from  a  hydrostatic  condition,  while  the  relationship  existing  between  the  average  stress 
(negative  of  pressure)  and  hydrostatic  strain  (volumetric  dilatation)  is  an  EOS. 


The  approach  proposed  by  Segletes  (34),  which  was  implemented  as  the  transversely  isotropic 
model  in  CTH  (Taylor  [729]),  has  found  reflection  in  the  modeling  of  an  EOS  for  orthotropic 
materials  proposed  by  Anderson  et  al.  (36).  At  this  point,  the  approach  to  describe  an  EOS  for 
orthotropic  materials  proposed  by  Anderson  et  al.  (36)  is  discussed.  The  final  expression  for  the 
EOS  (for  the  general  orthotropic  case)  based  on  equations  5,  42,  and  55,  and  proposed  by 
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Anderson  et  al.  (36),  can  be  written  in  the  following  form: 


p  =  H- 


:t* 


+  pTe 


3j3 


3/3 


{  Ex  ( 1  VyZ  UZy  )  3~  Ex  (  V.yX  -\~  Vzx  Vyz  )  A  Ex  ( Vzx  +  Vyx  V zy  )  }  '  ^xx 
{  Ex  ( l^yx  A  VzxVyz  )  +  Ey(l  ^XZ  ^zx  )  +  Ey(uzy  +  ^xy^zx)}  *  £yy 
{Ex(  V zx  “I-  VyXVZy  )  E  Ey(  V zy  H“  L'xyL'zx  )  +  Ez{  1  l',xyl',yx)  }  '  ^zz 


(59) 


where  P^  is  the  anisotropic  Hugoniot  pressure.  The  first  important  difference  between 
Anderson’s  and  the  conventional  isotropic  approach  is  in  the  way  the  Hugoniot  pressure  is 
approximated.  The  conventional  approximation  of  the  Hugoniot  pressure  Ph  may  assume  cubic 
least  squares  curve  fit  in  p  as 


Ph  = 


Ai/i  +  A2p~  +  A-3p^  ,  p  >  0 
Aiji  ,  p  <  0 


(60) 


where  parameters  At,  i  =  1,  2,  3  are  determined  by  fit  to  experimental  shock  compression 
(Hugoniot)  data.  Note  that  for  an  isotropic  material,  A1  is  the  bulk  modulus.  Anderson  et 
al.  (36)  proposed  that  in  order  to  have  consistency  and  correct  stresses  in  the  elastic  regime  for  an 
orthotropic  material,  the  Hugoniot  pressure  P^  should  be  approximated  as 


PA  — 


A) /i  +  A2p~  +  AsjP  ,  p  >  0 


ri'ip  , 


p  <  0 


(61) 


where  Al  is  defined  by  the  following  quantity  (Lomakin  [728];  Anderson  et  al.  [id]): 


A  |  (Ex(  1  VyZUZy)  T  Ey /(I  Vxz^zx)  +  Ez(  1  VxyVyx) 

T  2  •  [ Ex(uyx  T  azxiSyZ)  T  Ex ( vzx  T  VyXuZy)  T  Ey(iszy  T  a zx a Xy )  ] ) 


(62) 


where  (3  =  1  —  vxyvyx  —  vxzvzx  —  vyzvzy  —  2vyxvzyvxz.  Anderson  et  al.  (36)  interpreted  A\  as  an 
“effective”  or  “average”  bulk  modulus  Ka  (first  effective  bulk  modulus  equation  56).  It  is 
important  to  note  that  expression  for  the  EOS  can  likewise  be  obtained  using  equation  57, 
interpreting  A\  as  the  second  effective  bulk  modulus.  The  A2,  A3  parameters  are  determined 
from  the  fitting  of  the  experimental  data  in  all  cases.  This  provides  an  appropriate  description  of 
material  behavior  at  high  pressures  and  reduces  to  the  correct  relations  at  small  volumetric  strains. 
Alternatively,  A\,  A2,  A3  can  be  analytically  determined  through  a  Taylor’s  series  expansion  of 
the  Hugoniot  pressure  P.  Assume  that  the  linear  approximation  between  the  shock  velocity  Us 
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and  particle  velocity  up  (equation  6,  52  =  53  =  0)  exists,  and  that  the  Us~up  intercept  c  equals  the 
ambient  bulk  speed  of  sound  c0.  It  follows  that  the  Hugoniot  curve  and  Taylor’s  series  expansion 
of  the  Hugoniot  pressure  P  with  respect  to  //  can  be  written  in  the  form 


A  Po  ■  cl  ■  n  •  (1  +  fl) 

H  [l-p-(S-l)}2  ’ 

A  =  PoCq 

A2  =  A!x  ■  [1  +  2(5  —  1)] 

A3  =  A[-  [2(5-1)  +  3(5-  l)2] 


U  =  cq  +  Sup 

where  c0,  the  ambient  bulk  speed  of  sound,  has  the  following  definition 


(63) 


(64) 


(65) 


Some  attempts  to  use  Anderson’s  model  to  simulate  the  behavior  of  composite  materials  under 
shock  loading  have  been  made  by  Chen  et  al.  (130),  Hayhurst  et  al.  (131-133),  and  Hiermaier  et 
al.  (134).  The  work  of  Hayhurst  et  al.  (131-133)  was  directed  towards  numerical  model 
development  for  the  Nextel  and  Kevlar/Epoxy  materials  subject  to  hypervelocity  impact.  They 
also  performed  the  experimental  inverse  flyer  test  (IFPT)  for  Nextel  and  Kevlar/Epoxy.  Their 
models  were  to  be  macro-mechanically  based  and  suitable  for  implementation  into  a  hydrocode 
coupled  with  EOS.  Each  layer/weave  of  the  material  was  not  to  be  modeled  explicitly  but 
represented  by  an  equivalent  volume  with  properties  on  a  macroscale,  which  are  representative  of 
the  combined  micro-mechanical  response  of  the  volume  of  the  material  under  the  loading 
conditions  considered.  The  Kevlar/Epoxy  IFPT  experimental  data  and  simulation  results  are 
presented  in  figure  4.  A  polynomial  EOS,  equations  59  and  60,  was  used  in  conjunction  with 
damageable  orthotropic  stiffness.  The  samples  recovered  from  impact  tests  gave  evidence  of 
phase  changes  during  the  impact  events  (figure  4).  The  Kevlar/Epoxy  samples  recovered  after  an 
impact  at  572  m/s  showed  a  reduced  thickness  with  a  residual  bending  strength.  Separate  Kevlar 
sheets  were  found  in  the  impact  vessel  during  impact  test  at  788  m/s.  These  Kevlar  sheets 
showed  a  vanishing  bending  strength  with  all  Epoxy  impregnation  to  be  evaporated.  Impact  test 
at  1015  m/s  produced  a  fine  Kevlar  dust  all  over  the  impact  vessel.  At  this  impact  velocity,  Kevlar 
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Figure  4.  The  Kevlar/Epoxy  IFPT  simulated  and  experimental  back  surface  velocities  for  572,  788, 
and  1015  m/s.  The  experimental  data  Kevlar/Epoxy  materials  recovered  after  flyer  plate 
testing  were  taken  from  C.  J.  Flayhurst,  Multi-Physics  Analysis  of  Hypervelocity  Impact: 
Successes  and  Challenges,  FENET  Presentation  (Noordwjik,  2003)  (133). 


also  undergoes  a  phase  change  of  thermal  decomposition.  The  weave  and  yarns  decompose  into 
a  dark,  wool-like  material.  It  is  important  to  note  that  the  average  slope  of  the  initial  rise  in  back 
surface  velocity  is  well  represented  along  with  the  sharp  rise  in  velocity  (figure  4)  for  all  three 
impact  velocities.  The  observed  differences  between  the  simulations  and  experiments  can  be 
attributed  partly  to  the  selected  EOS  approach  and  corresponding  damage  model  (i.e.,  compaction 
behavior  and  initial  density).  However,  Hayhurst  et  al.  (131,  132 )  models  provide  a  significant 
improvement  over  the  simulation  with  the  standard  orthotropic  material  models. 

It  was  mentioned  previously  that  in  the  case  of  an  isotropic  material,  the  hydrostatic  stress  (or 
pressure)  induces  a  change  of  scale,  while  the  deviatoric  stress  only  induces  a  change  of  shape.  It 
is  obvious  that  in  order  to  keep  this  property  for  anisotropic  materials  (e.g.,  for  orthotropic 
materials),  the  definition  of  “pressure”  needs  to  be  generalized,  because  a  hydrostatic  pressure 
(isotropic  state  of  stress)  applied  to  an  anisotropic  material  results  in  an  anisotropic  state  of  strain. 
In  other  words,  this  loading  will  result  not  only  in  a  change  of  scale,  but  also  in  change  of  shape. 
This  is  inconsistent  with  the  definition  of  the  “generalized  pressure.” 
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3.5  Generalized  Decomposition  of  Anisotropic  Constitutive  Relationship  Suitable  for 
Shock- Wave  Modeling 

The  definition  of  pressure  equation  2  leads  to  the  invariant  quantity  (the  contraction  of  the  stress 
tensor  and  the  unit  tensor  is  divided  by  the  norm  of  the  unit  tensor).  This  is  not  the  case  when  the 
pressure  is  defined  as  the  portion  of  the  mean  stress  that  varies  directly  with  the  volumetric  strain, 
which  is  not  an  invariant  quantity  (i.e.,  this  quantity  cannot  be  expressed  as  contraction  of  the 
stress  tensor  and  another  second  order  tensor).  However,  the  high-pressure  thermodynamic  EOS 
of  anisotropic  materials  must  be  modified  to  account  correctly  for  the  elastic  behavior  at  small 
volumetric  strain.  The  generalized  decomposition  of  the  stress  tensor  was  constructed  by 
Lukyanov  (37-41,  135,  136),  assumed  and  given  by  the  following  theorem: 

Theorem:  For  anisotropic  materials  of  any  symmetry,  where  stress  state  and  strain  state  are 
coupled  via  the  generalized  Hooke’s  law  al:)  =  =  Jijki&ij  subject  to  the  constraints 

for  elements  of  stiffness  and  compliance  matrixes  ( written  in  Voigt  notation )  such  as 

3  3 

(Cfc  1  +  Ck2  +  Ckz)2  0  ,  (Jkl  +  Jk2  +  Jk3t)2  7^  0  j  (66) 

k= 1  k= 1 

there  is  one  and  only  one  state  of  stress,  ffj  =  patp  that  results  in  only  volumetric  deformation, 

6ij  =  evSij,  as  defined  by  the  following  decomposition: 


&ij  P  Oiij  T-  S ij  ,  p  QtijSij  0  ,  (XijQLij  3  ,  O } J  0  V  l  f  j 


(67) 


p*  =  p  +  p" 


P  = 


PijO'ij 


pS  _  PijSij 


cj  ij  pai3 


(68) 


fiklOlkl  Pkl&kl 

where  alJ  and  /3ij  are  the  first  and  second  generalizations  of  the  Kronecker  delta  symbol,  p  is  the 
pressure  related  to  the  volumetric  deformation,  ps  is  the  pressure  related  to  the  generalized 
deviatoric  stress,  p*  is  the  total  generalized  pressure,  Sij  is  the  generalized  deviatoric  stress 
tensor,  ev  =  eijSij  is  the  volumetric  deformation,  Cijki  is  the  stiffness  matrix,  Jijki  is  the 
compliance  matrix,  Cij  =  CJt  are  elements  of  the  stiffness  matrix  ( written  in  Voigt  notation),  and 
Jij  =  Jji  are  elements  of  the  compliance  matrix  ( written  in  Voigt  notation). 


It  is  important  to  note  that  some  attempts  to  construct  generalized  decomposition  of  stress  tensor 
were  made  by  Lomakin  (128)  and  Sawyer  (137).  Lor  anisotropic  materials,  the  total  hydrostatic 
“pressure”  has  been  defined  {37-41,  135,  136)  and  given  as 


where  p  is  the  pressure  related  to  the  volumetric  deformation  and  S,3  is  the  generalized  deviatoric 
part  of  the  stress  tensor.  The  relation  equation  69  is  the  correct  generalized  “pressure”  for  the 
elastic  regime.  To  provide  an  appropriate  description  of  behavior  for  general  anisotropic 
materials  at  high  pressures,  an  EOS  for  p  (pressure  related  to  the  volumetric  deformation)  has  to 
be  defined.  The  Mie-Griineisen  EOS  equations  8  and  9  was  used  for  pEOS  by 
Lukyanov  (37-41,  135,  136)  to  describe  the  thermodynamic  (EOS)  material  response,  which  also 
describes  correctly  the  material’s  behavior  at  small  volumetric  strains.  Therefore,  an  appropriate 
description  of  general  hydrostatic  “pressure”  at  high  pressures  has  the  following  form: 

p*  =  pEOS  +  &&  .  (70) 

C%j  Pij 

Note  that  the  methodology  described  above  can  be  applied  for  all  anisotropic  materials  and 
represents  a  mathematically  consistent  generalization  of  the  conventional  isotropic  case.  The 
methodology  for  calculation  of  components  of  the  tensor  ai3  has  been  previously 
defined  (37-41,  135,  136).  The  elements  of  tensor  a,3  can  be  written  in  the  following  form: 

c^n  =  (Gn  +  C\2  +  C13)  •  3 Kc 


CC22  —  {C12  +  C22  +  C23)  •  3 Kc  (71) 

a33  —  (C*i3  +  C'23  +  C33)  •  3 Kq 


Kc  —  1/ \J 3  •  [(Cli  +  C\2  +  G13)  +  (C'12  +  C22  +  C23)  +  (C*13  +  C*23  +  G33)" 


Kc  = 


9  Kc 

_ II  1 1 2 _  o 

OiijOiij  =||  a  ||  =  o  , 


(72) 


(73) 


where  C%3  is  the  stiffness  matrix  (written  in  Voigt  notation).  The  relation  equation  73  describes 
the  norm  of  tensor  al3,  which  reduces  to  Kroncker’s  delta  symbol  norm  in  the  limit  of  isotropy. 
Therefore,  the  norm  of  ai3  is  taken  to  be  \/3.  Furthermore,  the  following  set  of  equations  for  the 
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components  of  the  tensor  Aj  can  be  written  (Lukyanov  [37-41,  135,  136]): 

Ali  =  ( J\i  +  4 i 2  +  A3)  •  3A s 
$22  =  {J\2  +  J22  +  A3)  •  3AS 


A3  —  (4 \:i  +  A3  +  A3)  '  3 Ks 


(74) 


As  —  1/ ^3  •  [(Ai  +  A2  +  A3)2  +  (A2  +  4 22  +  A3)2  +  (As  +  A3  +  J33Y 


PijPij  =11/5  II 2=  3 


where  4tJ  are  elements  of  compliance  matrix  (written  in  Voigt  notation).  To  be  consistent  with 
the  definition  of  the  bulk  speed  of  sound  c  for  isotropic  material,  the  following  definition 


c  = 


(75) 


is  assumed  for  anisotropic  bulk  speed  of  sound.  Here  the  first  generalized  bulk  modulus  Kc  is 
defined  according  to  equation  72.  Based  on  methodology  previously  described,  we  can  conclude 
that  the  two  fundamental  tensors  and  which  represent  material  properties,  have  been 
defined.  Both  of  them  can  be  considered  as  generalizations  of  the  Kronecker  delta  symbol,  which 
plays  the  main  role  in  the  theory  of  isotropic  materials.  Using  two  fundamental  tensors  atJ  and 
Aj,  the  definitions  of  “total  pressure”  and  pressure  corresponding  to  the  volumetric  deformation 
can  be  defined.  In  the  limit  of  isotropy,  tensors  al3  and Aj  have  the  following  values 
an  =  «22  =  «33  =  1,  Ai  =  P22  =  A 3  =  1  and  the  proposed  generalization  returns  to  the 
traditional  classical  case,  where  tensors  atj  and  Ai  equal  d%]  and  equations  67,  68,  and  69  take  the 
following  form: 


O' kk  4ijO~ij 

Q%]  3  Pij&ij 


&kk 

IT 


P*  =P  ,  Sij 


Sij  = 


a , 


O'kk 

IT 


•  (76) 


Here,  p*  =  p  is  the  conventional  hydrostatic  pressure  and  StJ  is  the  conventional  deviatoric  stress 
tensor.  Also,  the  two  parameters  Kc  and  Ks  were  considered  as  the  first  and  second  generalized 
bulk  moduli.  In  the  limit  of  isotropy  they  reduce  to  the  well-known  expression  for  the 
conventional  bulk  modulus  K  =  E /  [3(1  —  2u)\. 


Most  of  the  work  discussing  the  response  of  composites  to  shock  loading  has  examined  the 
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material  with  the  fiber  plane  normal  to  the  loading  axis  ( e.g .,  Dandekar  et  al.  [13],  Zhuk  et 
al.  [73S],  Riedel  et  al.  [139],  Zaretsky  et  al.  [140]).  Overall,  when  shocked  in  this  orientation, 
fiber-based  composites  appear  to  behave  in  a  manner  similar  to  monolithic  polymers  (e.g.,  Carter 
et  al.  [141],  Millett  et  al.  [142]),  indeed  Zaretsky  et  al.  (140)  have  made  this  point.  The 
experimental  study  of  a  carbon  fiber  epoxy  composite,  shocking  along  the  through-thickness 
orientation  axis,  showed  no  evidence  of  an  inelastic  deformation.  Therefore,  using  an  anisotropic 
nonlinear  continuum  framework  and  generalized  decomposition  of  a  stress  tensor 
equations  66-73,  shock-wave  propagation  in  CFC  materials  was  examined  (Lukyanov  [41]).  The 
results  of  comparison  between  experimental  data  and  numerical  simulation  are  shown  in  the 
figure  5.  The  important  characteristic,  the  arrival  time  to  the  Hugoniot  stress  level  at  the  0  mm 
position  and  back  surface  are  in  good  correlation  with  experimental  data.  Further  comparison 
shows  that  the  pulse  width  and  the  reloading  trace  are  in  good  agreement  with  the  experimental 
data.  The  maximum  difference  between  the  experimental  data  and  new  proposed  model  for  the 
plateau  stress  was  6%. 


Time  ((is)  Time  (ps) 

Figure  5.  Representative  experimental  gauge  traces  from  the  through-thickness  orientation  at  the 

0  mm  position  and  at  the  back  surface,  respectively  (see  Millett  et  al.)(15).  The  specimen 
was  3.8  mm  thick.  The  impact  conditions  were  a  5  mm  dural  flyer  at  V  =  504  m/s.  The 
dotted  curve  is  the  numerical  data  obtained  using  proposed  damage  model;  the  solid  curve  is 
the  experimental  data. 

These  experimental  results  show  that  the  relationship  between  shock  velocity  and  particle  velocity 
through  the  thickness  orientation  is  linear,  yielding  the  relation  (Millett  et  al.  [75]): 
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Ug  =  Cq  +  S[‘up,  where  Cp  =  3230  m/s,  S[  =  0.92.  The  experimental  data  in  figure  6  are 
obtained  for  target  plate  thicknesses  in  the  range  between  2.3  and  5.7  mm  (see  Millett  et  al.  [75]). 
It  can  be  seen  that  there  is  a  degree  of  scatter  in  this  data.  It  was  noticed  that  EOS  data  Cq  and 
S[J  are  insensitive  to  thickness  of  the  target  plate  over  experimental  range  of  the  investigation  (see 
Millett  et  al.  [75]). 


0  0.2  0.4  0.6  0.8  1 

Particle  Velocity  (mm  jjs'1) 


Figure  6.  Experimental  data  Ug-up  for  the  carbon- liber-composite  material,  showing  the  variation  with 
specimen  thickness  (experimental  data  obtained  by  Millett  et  al.)(15).  The  dotted  curve  is 
calculated  using  experimental  data  for  Ug  -up\  the  solid  curve  is  calculated  using  numerical 
simulation  based  on  the  material  model  equations  66-73. 

However,  it  would  be  expected  that  when  the  fiber  direction  is  orientated  along  the  impact  axis, 
the  response  would  be  quite  different.  The  results  of  Eden  et  al.  (90)  clearly  differentiated 
between  the  response  of  the  fibers  and  the  matrix,  indicating  that  the  individual  fibers  were  acting 
as  wave-guides.  Bordzilovsky  et  al.  (12)  examined  the  effects  of  orientation,  with  the  orientation 
of  the  fibers  ranging  from  5°  to  90°  to  the  shock  axis.  Where  the  mis -orientation  between  fibers 
and  the  loading  axis  was  small,  a  distinct  low  amplitude  precursor  wave  was  observed  before 
arrival  of  the  main  shock.  As  the  angle  increased,  the  duration  of  this  precursor  decreased  until  it 
disappeared  at  90°.  This  was  interpreted  as  an  elastic  wave.  Hereil  et  al.  (11)  observed  similar 
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behavior  in  a  three-dimensional  carbon-carbon  composite,  with  a  low  amplitude  ramp  preceding  a 
much  more  quickly  rising  shock  front.  However,  their  interpretation  suggested  that  this  precursor 
was  due  to  a  high  velocity  wave  transmitted  along  the  fibers  orientated  in  the  shock  axis,  while  the 
main  shock  was  transmitted  through  the  “matrix”  (figure  7).  Millett  et  al.  (15)  examined  the 
carbon  fiber-epoxy  composite  and  obtained  that  the  trace  clearly  has  a  ramped  nature  (figure  7),  in 
common  with  the  observations  of  Eden  et  al.  (90)  when  using  quartz  stress  gauges  to  measure  the 
shape  of  the  stress  pulse  in  a  similar  orientation.  This  is  qualitatively  similar  to  the  traces 
observed  by  Hereil  et  al.  (11)  and  Bordzilovsky  et  al.  (12),  respectively. 


Figure  7.  Representative  experimental  gauge  traces  from  the  fiber  0°  orientation  (see  Millett  et  al.)(15). 
The  specimen  was  10  mm  thick.  The  impact  conditions  were  a  5  mm  copper  flyer  at 
V  =  936  m/s.  According  to  Hereil  et  al.  (11).  the  precursor  was  due  to  a  high  velocity  wave 
transmitted  along  the  fibers  orientated  in  the  shock  axis,  while  the  main  shock  was  transmitted 
through  the  matrix. 


It  was  pointed  out  by  Bethe  (143)  that,  for  stable  shock  waves,  the  shock  velocity  must  increase 
with  pressure.  Tins  means  that  if  the  shock  velocity  should  decrease  with  pressure,  then  the 
shock  front  would  break  up  into  two  or  more  waves,  or  possibly  one  wave  with  a  continuously 
smeared  front.  Note  that  the  experimental  data  for  the  through  thickness  orientation  (see, 
Lukyanov  [775])  can  be  fitted  by  a  linear  relation,  and  there  is  no  explicit  evidence  of  the  shock 
front  breaking  up;  however,  the  analysis  of  the  experimental  data  of  Millett  et  al.  (15)  for  selected 
CFC  material  shows  that  the  shock  velocity  along  the  fiber  0°  orientation  decreases  with 
pressure — therefore,  a  two-wave  structure  is  proposed  for  describing  the  experimental  data. 
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Furthermore,  as  the  severity  of  the  shock  increases,  the  Hugoniot  stress  levels  of  the  two 
orientations  converge.  This  fact  demonstrates  that  the  selected  CFC  material  shows  isotropic 
behavior  at  high  shock  intensities,  and  can  be  described  as  an  isotropic  mixture  of  epoxy  binder 
and  fractured  fibers  (see,  Lukyanov  [755]).  This  is  qualitatively  similar  to  the  Kevlar-Epoxy 
behavior  observed  by  Hayhurst  et  al.  (131-133)  and  Hiermaier  el  al.  (134)  (figure  4). 


4.  Anisotropic  Plasticity  Flow 


The  anisotropic  material  response  under  shock  loading  often  shows  inelastic  (e.g.,  plastic) 
deformation.  The  ability  to  predict  inelastic  response  can  be  important  for  both  ductile  metals 
and  relatively  brittle  composite  materials.  Note  that  even  though  composite  materials  are 
generally  considered  to  be  brittle,  it  is  often  prudent  to  use  a  small  value  for  the  equivalent  plastic 
strain  as  a  failure  criterion  (i.e.,  that  point  at  which  the  material  can  no  longer  support  shear  and 
tensile  stresses).  Traditionally,  the  incremental  strain  of  a  material  element  is  written  as  the  sum 
of  elastic  and  inelastic  (e.g.,  plastic)  strains  AetJ  =  Aefj  +  Ae*”.  Two  types  of  inelastic 
deformation  models,  which  do  not  result  in  loss  of  cohesion,  are  considered  in  the  shock-wave 
modeling  in  anisotropic  materials:  an  anisotropic-yield-criteria-based  plasticity  model 
(associated  and  non-associated)  and  a  dislocation-based  plasticity  model. 

4.1  Associated  and  Non-associated  Anisotropic  Plasticity 

The  main  aspects  of  a  phenomenological  constitutive  model  can  be  characterized  by  a  yield 
criterion  representing  a  surface  that  separates  the  elastic  and  plastic  regions  of  the  stress  space,  a 
flow  potential  gradient  that  represents  the  direction  of  plastic  strain  rate,  and  a  strain  hardening 
rule.  Generally,  a  yield  function,  a  flow  rule,  and  a  set  of  evolution  equations  for  M  state 
variables  are  required  to  fully  describe  the  constitutive  relations  of  a  plastic  material.  The  yield 
function  F  for  an  anisotropic  material  in  classical  formulation  can  be  expressed  in  the  following 
general  form: 

=  0,  k  =  1,  ...,M  .  (77) 

Here  at]  are  the  stress  components  and  represents  a  set  of  state  variables.  The  subscript  k  is 
introduced  to  indicate  that  there  may  be  several  state  variables  including  the  hardening 
parameters.  When  the  material  deforms  plastically,  the  inelastic  part  of  deformation  (plastic 
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deformation)  is  defined  by  the  flow  rule 


de?-  =  d\-^~  ,  (78) 

O  (7  ij 

where  g  =  g(ctij,  7fc)  is  the  plastic  potential,  del  are  the  differentials  of  the  plastic  strain 
components,  and  dX  is  a  positive  scalar.  In  this  study,  we  consider  the  associated  plasticity  flow 
models  ( i.e .,  F  =  g)  and  non-associated  plasticity  flow  models  (i.e.,  F  ^  g).  Finally,  the 
evolution  of  M  state  variables  can  be  described  by 

oik  =  Ak (de^j ,ai:j,g/i)  ,  k  =  1, M;  l  =  1, M  .  (79) 

For  complex  plasticity  models,  several  evolution  equations  may  be  defined  and  the  forms  of  these 
equations  can  be  very  complex. 

Several  anisotropic  yield  criteria  and  their  associated  yield  surfaces  have  been  developed  in  the 
past.  These  criteria  have  included  a  maximum  stress  law,  a  maximum  strain  law,  and  quadratic 
laws.  One  of  the  first  yield  conditions  proposed  for  the  anisotropic  material  is  the  quadratic  yield 
criterions  proposed  by  Dorn  (144)  and  von  Mises  (Hill  [145])  for  plastically  incompressible  metal 
crystals  with  different  lattice  symmetry.  In  the  theory  of  anisotropic  yield  criteria,  the  most 
well-known  work  is  HilFs  quadratic  formulation  (Hill  [145]),  which  contains  six  parameters 
specifying  the  state  of  anisotropy,  but  is  similar  in  form  to  the  Mises ’s  criterion  for  isotropic 
metals.  From  the  literature  review,  a  general  group  of  anisotropic  yield  criteria  suitable  for 
anisotropic  metals  can  be  found.  This  group  includes  the  yield  conditions  proposed  by 
Bassani(i46)  and  Hosford(i47);  the  yield  surface  of  four  degrees  specified  by  Gotoh(748)  and 
Arminjon  et  al.  (149);  as  well  as  yield  surfaces  of  k  degrees  analyzed  by  Barlat  and  Lian  (150). 
Barlat  et  al.  (151,  152),  Karafillis  and  Boyce  (153),  Barlat  el  al.  (154,  155),  Bron  and 
Besson  (156),  Darrieulat  and  Montheillet  (157),  Stoughton  and  Yoon  (758),  Kowalczyk  and 
Gambin  (159),  Hu  (160),  Hashiguchi  (161),  Hu  (162,  163),  and  Barlat  et  al.  (164).  The  quadratic 
laws  are  robust  and  are  particularly  well  suited  to  multiaxial  stress  state.  In  particular,  the 
Tsai-Hill  theory  (Hill  [145];  Tsai  and  Hahn  [i<55])  was  considered  by  Segletes  (34),  Anderson  et 
al.  (36),  and  De  Vuyst  et  al.  (166). 

Several  models  for  inelastic  deformation  have  been  proposed  based  on  generalized  decomposition 
of  stress  tensor  by  Lomakin  (128)  and  Sawyer  (137).  The  mathematically  consistent  yield 
function  of  a  fully  anisotropic  material  based  on  generalized  decomposition  of  the  stress  tensor 
into  generalized  spherical  part  (generalized  hydrostatic  pressure)  and  generalized  deviatoric  stress 
tensor  equation  67  was  proposed  by  Lukyanov  (37-41).  Based  on  research  of  Spitzig  and 
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Richmond  (167),  Stoughton  and  Yoon(75S)  and  generalized  decomposition  equation  67,  the 
following  yield  function  was  used  for  modeling  anisotropic  elastic-plastic  shock  waves: 

F(Sij)  =  1  +  XP*)  <  Y(ep)  ,  f  =  <J‘ja'J  =  yA T2^ijotij  ,  P*  =  ~V*  ,  (80) 

OLklOLkl  ||  OL  || 

where  p*  is  the  generalized  hydrostatic  stress  and  T  ( S,3 )  is  described  by  generalized  Hill’s  yield 
function: 


—  F(a33Syy  —  a22Szz)2  +  G(a>uSzz  —  ai33Sxx )2  +  H(a22Sxx  —  an  Syy)2 
+  2  NS%  +  2  LSI  + 

The  material  constants^;,  F,  G,  H,  N,  L,  and  M  are  specified  in  terms  of  selected  initial  yield 
stresses  in  uniaxial  tension,  compression,  and  equibiaxial  tension.  It  is  worth  noting  that 
plasticity  model  equation  78  is  naturally  independent  from  the  generalized  hydrostatic  stress 
p*  =  ( CTjjCCy ) / (aijaij),  and  therefore,  the  following  equality  can  be  written: 

^2(Sij)  =  ^2(aij)  ,  o,j  pP\jj  -  St;  ,  p*  ^  0  ,  (82) 

where  Sl3  is  the  generalized  deviatoric  stress  tensor  and  al3  is  the  stress  tensor.  This  yield 
function  was  validated  for  a  number  of  materials,  e.g.,  AA2008  T4,  AA2090  T3,  AA7108  Tl,  and 
AA6063  Tl  (Lukyanov  [40]).  The  shock-wave  propagation  in  the  anisotropic  aluminum  alloy 
7010-T6  using  the  yield  function  equation  78  was  performed  by  Lukyanov  (39). 

The  experimental  values,  0.39  and  0.33  GPa,  for  elastic  response  from  the  longitudinal  and  short 
transverse  directions,  respectively  (figures  8  and  9),  are  in  good  correlation  with  the  modeled 
values  of  the  HEL  longitudinal,  0.395  GPa,  and  short  transverse,  0.333  GPa.  The  errors  with 
respect  to  the  experimental  values  are  1.4%  and  0.9%,  respectively,  to  the  longitudinal  and  short 
transverse  directions.  Besides,  other  important  characteristics,  the  arrival  time  to  the  HEL  and 
the  plastic  wave  velocity,  are  in  good  correlation  with  experimental  data.  Further  comparison 
shows  that  the  pulse  width  and  the  reloading  trace  are  in  good  agreement  with  the  experimental 
data  (figures  8  and  9).  The  main  conclusion  obtained  from  these  results  is  that  the  non-associated 
anisotropic  plasticity  model,  as  it  stands,  is  suitable  for  simulating  elasto-plastic  wave 
propagation  in  anisotropic  solids.  Besides,  different  HELs  are  obtained  when  the  material  is 
impacted  in  different  directions;  their  excellent  agreement  with  the  experiment  demonstrates  that 
the  anisotropic  plasticity  model  is  adequate. 

However,  further  work  is  required  both  in  the  experimental  and  constitutive  modeling  areas  to 
find  a  better  description  of  anisotropic  material  behavior. 
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Plate  Impact  Test  —  Longitudinal  Direction 


8.  Back-surface  gauge  stress  traces  (longitudinal  direction)  from  plate-impact 
experiments  vs.  numerical  simulation  of  stress  (PMMA)  waves  for  plate  impact 
tests  (impact  velocities  450  and  895  m/s)  -  target  AA7010  T6. 


Plate  Impact  Test  —  Transverse  Direction 


Figure  9.  Back-surface  gauge  stress  traces  (transverse  direction)  from  plate-impact 

experiments  vs.  numerical  simulation  of  stress  (PMMA)  waves  for  plate  impact 
tests  (impact  velocities  450  and  895  m/s)  -  target  AA7010  T6. 
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4.2  Dislocation-based  Plasticity 


Much  of  the  dislocation  theory  is  now  well  established  and  is  attributed  to  the  early  efforts  of 
Taylor!/ 6#),  Orowan (769),  Gillis  and  Gilman  (170,  171),  and  Gilman  (172,  173).  To  incorporate 
dislocation-based  plasticity,  the  Orowan  equation  (Gilman  [775];  Winey  and  Gupta [28,  29])  is 
used  to  relate  the  motion  of  dislocations  on  a  slip  plane  (denoted  by  a)  to  the  plastic  shear  strain 
rate  on  that  plane: 


la  m  r 


(83) 


where  iV,“  is  the  mobile  dislocation  density  operative  on  the  slip  plane,  ba  is  the  magnitude  of  the 
Burgers  vector,  and  Va  is  the  average  dislocation  velocity.  The  incremental  plastic  strains  in  the 
coordinate  system  used  to  describe  wave  propagation  are  calculated  by  appropriately 
transforming  and  summing  the  plastic  shear  strains  from  all  the  individual  slip  systems  (Winey 
and  Gupta  [28-30]): 

=  \  (a“*a3i  +  °3 i°lj)  7a At  ,  (84) 

G! 

where  At  is  the  time  step  encountered  in  wave  propagation  calculation,  and  a"  is  the  rotation 
matrix  relating  the  coordinate  system  used  to  describe  wave  propagation  to  the  slip  plane 
coordinates.  The  resolved  shear  stress  ra  causing  dislocation  motion  for  a  given  slip  system  is 
determined  by  transforming  the  Cauchy  stress  from  the  coordinate  system  used  to  describe  wave 
propagation  to  the  slip  plane  coordinates  system: 


^a  _  _a  _  „a  _ 

r  —  cr13  —  aua3jaij 


(85) 


To  account  for  the  rotation  of  slip  planes  due  to  finite  elastic  strains,  the  rotation  matrices  a"  are 
updated  each  moment  of  time  through  a  transformation  that  depends  on  the  incremental  elastic 
deformations  (Winey  and  Gupta  [28-30]).  The  dislocation-based  models  produce  only  shear 
strains,  resulting  in  plastic  incompressibility  Depkk  =  0.  The  mobile  dislocation  density  TV"  and 
the  average  dislocation  velocity  Va  must  be  defined  in  order  to  perform  numerical  simulations 
using  equation  77.  Winey  and  Gupta  (28-30)  employed  a  model  (Gilman  [173])  in  which 
dislocations  undergo  regenerative  multiplication  under  shock  loading  via  a  multiple  cross  glide 
mechanism.  With  this  mechanism,  the  dislocation  density  is  linearly  related  to  the  accumulated 
plastic  shear  strain  yj  (assumed  positive) 

N«  =  Nm0  +  My£  ,  (86) 

where  Nm0  is  the  initial  density  of  mobile  dislocations  and  M  is  the  multiplication  parameter. 
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The  model  also  incorporates  a  stress-dependent  dislocation  velocity 


V°  =  Vo  exp  [—79/  (ra  —  r0)]  ,  (87) 

where  79  is  the  drag  stress  parameter,  rQ  is  the  resolved  shear  stress,  r0  is  the  threshold  stress  for 
dislocation  motion,  and  V0  is  the  shear  wave  speed.  The  dislocation-based  plasticity  model  was 
used  by  Winey  and  Gupta  (29)  for  a  6061-T6  aluminum  plate  impacting  the  LiF  crystals  at  a 
velocity  of  ca.  0.34  km/s,  where  the  target  crystals  are  backed  by  a  quartz  stress  gauge. 


5.  Early  Anisotropic  Implementations  in  Hydrocodes 


Computational  methods,  as  proposed  by  von  Neumann  (Clippinger  [174])  came  to  the  forefront 
with  the  advent  of  computers,  beginning  with  the  ENIAC,  designed  and  built  by  the  University  of 
Pennsylvania  (dedicated  1946)  and  installed  at  the  United  States  Ballistic  Research  Laboratory  in 
1947  for  the  calculation  of  projectile  trajectories  (Bergin  [775]).  By  the  mid  1950s,  newer 
generations  of  computers  were  fully  engaged  in  computational  fluid  dynamics,  e.g.,  the 
particle-in-cell  method  (Evans  and  Harlow  [176]). 

By  the  early  1960s,  the  hydrocodes  had  evolved  to  model  elasticity  and  radial-return  (isotropic) 
plasticity  (Wilkins  [102]).  It  was  not  until  the  1970s,  though,  that  elastic  and  plastic 
transverse-isotropy  were  computationally  implemented  in  the  hydrocode  framework  of  the  HELP 
code  (Sedgwick  et  al.  [177]).  Their  formulation,  however,  cannot  be  termed  deviatoric.  The 
form  of  the  relation  used  in  HELP  is 

A Sij  =  Cijkl Aekl  -  K(AV/V)8ij  ,  (88) 

where  K  is  the  bulk  modulus,  which  presumably  can  be  made  dependent  on  dilatation  (and 
therefore  hydrostatic  stress).  In  this  way,  the  formulation  may  also  provide  the  flexibility  akin  to 
a  truly  deviatoric  formulation.  However,  equation  88  is  not  a  deviatoric  relation,  since  the 
deviatoric  stress  increment  is  not  related  to  deviatoric  strain  increment,  but  rather  is  expressed  in 
terms  of  the  total  strain  increment.  Similarly,  the  bulk  modulus  (as  opposed  to  the  effective  bulk 
moduli  derived  in  equations  56  and  58)  is  functionally  dependent  on  deviatoric  stress,  and  in  this 
sense,  equation  88  exhibits  flawed  behavior  if  the  deviatoric  variation  in  bulk  modulus  is  not 
modeled.  Finally,  the  flexibility  afforded  in  equation  88  by  allowing  the  bulk  modulus  to  vary 
with  hydrostatic  stress  has  the  disturbing  effect  that  the  resulting  sum  of  the  normal  stress 
deviators  is  not  generally  zero.  Thus,  the  use  of  the  term  “stress  deviators”  to  describe  the 


38 


left-hand  side  of  equation  88  does  not  even  seem  justified. 


The  EPIC  code  added  anisotropy  in  1980  (Johnson  et  al.  [178])  and  uses  a  form  similar  to 
equation  88  except  that  K  is  defined  in  such  a  way  as  to  force  the  sum  of  the  normal  stress 
deviators  to  zero.  This  ad  hoc  procedure  will  coincidentally  mimic  the  behavior  of  equation  49, 
though  the  formulation  is  in  error  during  the  subsequent  hydrostatic  stress  calculation  by  not 
accounting  for  the  deviatorically  induced  dilatation  e  (see  equation  46).  In  this  sense,  both  the 
algorithms  of  HELP  and  EPIC  introduced  errors  that  were  sensitive  to  the  size  of  the  hydrostatic 
strain  increment  (i.e.,  the  computational  timestep)  if  variable  compressibility  were  employed  in 
the  EOS.  This  problem  was  alleviated  by  the  algorithm  proposed  by  Segletes  in  1987  (34),  as 
described  in  section  3.4  of  this  report.  Segletes’  model  was  implemented  into  the  Lagrangian 
DEFEL  code  and  was  used  to  study  the  effects  of  shear  forming  of  copper  liners  upon  the  collapse 
of  shaped-charge  warheads  (179,  180). 

Some  early  computational  implementations  of  anisotropy  restricted  themselves  to  linear  elastic 
anisotropy,  thereby  limiting  their  application  to  dynamic  problems  involving  gentle  non-plastic 
impact.  For  example,  the  TIP  code  of  Zak  and  Pillasch  (181)  is  in  this  category,  as  well  as  the 
anisotropic  implementation  for  the  TOODY  code  by  Swegle  and  Hicks  (182).  Zak(183) 
introduced  anisotropic  plasticity  into  the  SANX  code,  thereby  allowing  calculation  of  a  wider 
class  of  impact  problems.  With  SANX,  while  the  plasticity  was  anisotropic,  the  elasticity 
appears  to  be  treated  in  an  isotropic  fashion.  Later,  the  orthotropic  elastic  model  was 
implemented  into  Lawrence  Livermore  National  Laboratory  (LLNL)  DYNA3D  and  NIKE3D 
public  domain  codes  (Hallquist  and  Whirley  [I84\:  Maker  et  al.  [785]),  and  in  LS-DYNA 
commercial  code  (Hallquist  [186]).  In  1995,  the  CTH  code  of  Sandia  National  Laboratories 
(Taylor  [129])  incorporated  the  transversely  isotropic  model  of  Segletes  (34)  as  its  “TI  model.” 

Nowadays,  some  models  incorporate  a  user-defined  subroutine  within  the  commercial  software 
(e.g.,  ABAQUS)  to  take  into  account  either  a  homogenous  orthotropic  model  that  examines  the 
bulk  macroscopic  deformation  response  or  the  discrete  non-homogeneous  material  distribution 
(Pankow  et  al.  [187]).  In  such  models,  EOS  does  not  even  play  a  role.  Rather,  the  volumetric 
treatment  is  one  of  a  constant  bulk  modulus,  insensitive  to  thermal  effects. 

Only  recently,  some  attempts  were  made  to  properly  account  for  nonlinear  shock  effects,  energy 
dependence,  and  anisotropic  stiffness  in  the  existing  hydrocodes.  The  approach,  based  on  that  by 
Anderson  (36),  has  been  incorporated  into  the  AUTODYN  hydrocode  (Hayhurst  [131]).  The 
generalized  decomposition  of  anisotropic  constitutive  relationship  and  EOS  suitable  for 
shock-wave  modeling  in  anisotropic  materials  have  been  incorporated  into  the  DYNA3D 
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hydrocode  (Lukyanov  [3§]). 


6.  Conclusion 


In  the  area  of  studies  of  materials  under  shock  loading,  barring  a  few,  most  of  the  studies  are  of  a 
continuum  nature.  Therefore,  the  questions  as  to  how  shocks  propagate  and  what  the  state  of 
matter  is  immediately  behind  the  shock  front  are  pretty  much  open  questions  for  anisotropic 
materials,  though  a  few  attempts  have  been  made.  It  is  expected  that  more  experiments  of 
shock-wave  propagation  will  provide  very  valuable  information.  The  plate  impact  tests  in 
different  directions  at  low  and  high  intensities  hold  significant  promise,  particularly  because  of 
their  one-dimensional  planar  shock-wave  propagation  and  their  time -resolved  nature.  Interaction 
of  theory  and  experiment  in  this  region  may  provide  insight  into  the  dynamic  stress  relation 
processes  behind  the  shock  front.  The  fact  that  the  shock  jump  in  pressure  is  essentially  a 
nonequilibrium  phenomenon  may  also  play  a  deciding  role.  In  the  context  of  general  shock- wave 
research,  the  measurement  of  precise  temperature  is  a  crying  need.  Presently,  even  the  improved 
pyrometric  method  gives  temperatures  that  are  substantially  higher  than  expected  even  at  low 
pressure  ranges.  The  area  of  study  of  kinetics  of  phase  transitions  of  anisotropic  materials  under 
shocks  is  an  interesting  and  challenging  one,  and  some  information  could  be  derived  from 
accurate  experimental  measurements.  On  the  theoretical  side,  first-principle  molecular  dynamic 
calculation  for  anisotropic  materials  may  contribute  significantly  to  the  understanding  of  shock 
propagation  at  high  pressures.  It  can  be  hoped  that  this  interaction  of  theory  and  experiments  will 
open  up  new  vistas  in  the  area  of  high-pressure  physics  of  anisotropic  materials.  Besides,  there  is 
an  interesting  possibility  that  the  strain-rate  sensitivity  is  itself  orientation  dependent  in 
anisotropic  materials.  It  has  been  noted  for  some  anisotropic  metals  that  the  spall  strengths  are 
similar  in  both  orientations  at  lower  impact  stresses,  while  at  higher  levels,  the  spall  strength  is 
higher  in  the  longitudinal  direction.  This  would  seem  to  indicate  a  higher  degree  of  strain-rate 
sensitivity  in  the  longitudinal  orientation,  and  would  seem  to  agree  with  the  observations  made 
with  the  HELs  for  some  anisotropic  materials.  Therefore,  further  development  of  the  constitutive 
equations  taking  into  account  strain  rate  sensitivity  is  required. 

Modem  hydrocode  shock  modeling  capabilities  are  confined  almost  exclusively  to  isotropic 
media;  little  provision  has  been  made  for  anisotropic  materials.  In  order  to  make  numerical 
simulations  of  a  hypervelocity  impact  on  the  reference  configuration  possible  in  terms  of 
computation  times,  a  macroscopic  continuum  model  for  the  involved  materials  had  to  be 
developed.  Due  to  the  experimentally  observed  behavior  of  anisotropic  materials,  new 
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orthotropic  hydrocode  models  were  developed  from  the  theoretical  basis.  This  was  necessary 
because  total  generalized  pressure  or  conventional  pressure  inside  these  materials  depends  on 
deviatoric  strain  components  as  well  as  volumetric  strain.  Nonlinear  effects,  such  as  shock 
effects,  can  be  incorporated  through  the  volumetric  straining  in  the  material.  Thus,  a  basis  was 
found  to  couple  the  anisotropic  material  stiffness  and  strength  with  an  anisotropic  shock  effects, 
associated  energy  dependence,  and  damage  softening  process.  To  our  knowledge,  this  report 
presents  the  current  state  of  the  art  in  the  experimental  and  theoretical  developments  of  this 
fascinating  field. 
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